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Abstract: Systems biology aims at achieving a system-level understanding of living 
organisms and applying this knowledge to various fields such as synthetic biology, 
metabolic engineering, and medicine. System-level understanding of living organisms can 
be derived from insight into: (i) system structure and the mechanism of biological networks 
such as gene regulation, protein interactions, signaling, and metabolic pathways; (ii) system 
dynamics of biological networks, which provides an understanding of stability, robustness, 
and transduction ability through system identification, and through system analysis methods; 
(iii) system control methods at different levels of biological networks, which provide an 
understanding of systematic mechanisms to robustly control system states, minimize 
malfunctions, and provide potential therapeutic targets in disease treatment; (iv) systematic 
design methods for the modification and construction of biological networks with desired 
behaviors, which provide system design principles and system simulations for synthetic 
biology designs and systems metabolic engineering. This review describes current 
developments in systems biology, systems synthetic biology, and systems metabolic 
engineering for engineering and biology researchers. We also discuss challenges and future 
prospects for systems biology and the concept of systems biology as an integrated platform 
for bioinformatics, systems synthetic biology, and systems metabolic engineering. 
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1. Introduction 

The Human Genome Project and high-throughput experimental methodologies such as microarray 
chromatin-immunoprecipitation DNA chips (ChlP-chip) have led to the development of biology as an 
increasingly information-rich science encompassing transcriptomes, proteomes, metabolomes, 
interactomes, and so forth [1,2]. Some have suggested that systems biology is nothing more than a new 
name for integrative physiology, which has been practiced for the past 50 years or more. Because of 
these novel technologies, biologists have been able to collect data at a rate that was unimaginable a 
decade ago. The context of biology has profoundly changed over the past 20 years. These changes 
provide a powerful new framework for systems biology that moves it far beyond classical integrative 
physiology. A systems biology approach implies that every system of any level of biological systems 
can be analyzed with respect to the system's structure, in particular, in terms of its dynamics, method of 
control, and method of system design. Systems biology involves genomic, transcriptomic, proteomic, 
and metabolic investigations from a systematic perspective. As a result, systems biology has become the 
frontier of modern biological research; large amounts of new omics data cannot be understood without a 
network or systems viewpoint and without highly sophisticated computational analyses [3-11]. 

The role of systems biology in modern biological research (Figure 1) requires powerful 
computational tools to mine large-scale data sets of information on genetics, proteins, DNA-protein 
binding, metabolism, and so forth. These tools are used to construct dynamic system models for the 
interpretation of specific mechanisms of some cellular phenotypse (or behaviors) from a system (or 
network) perspective [12-17]. To construct a dynamic system model of biological networks from omics 
data, system identification technologies (i.e., reverse-engineering schemes) are needed to estimate the 
parameter values of dynamic models and the order of biological networks [18-44]. Synthetic biology 
metabolic engineering has recently been developed for designing synthetic genetic networks for the 
production of specific cellular functions in host cells [45-51]. Based on system models and mechanisms 
in systems biology, synthetic genetic circuits and metabolic engineering pathways can be designed to 
investigate cellular behaviors [52-63]. These synthetic biological technologies can be employed to 
investigate the models and mechanisms of systems biology. Discrepancies between the real behavior of 
synthetic genetic networks and the desired behavior predicted by the models and mechanisms of systems 
biology can be fed back to modify the models through methodologies of systems synthetic biology and 
systems metabolic engineering. Based on the role of systems biology (Figure 1), this review describes 
current developments in bioinformatics, systems synthetic biology, and systems metabolic engineering. 
It discusses how systems biology can serve as an integrated platform for bioinformatics, systems 
synthetic biology, and systems metabolic engineering in the future. 

Bioinformatics is crucial in genome-wide analyses for understanding cell physiology at different 
cellular levels (i.e., genome, transcriptome, proteome, and metabolome levels) [12,13]. The various 
disciplines of bioinformatics provide invaluable information on the global cellular status for systems 
biology, systems synthetic biology, and systems metabolic engineering, as well as a thorough analysis of the 
cell. Genomic information in bioinformatics represents the whole genetic makeup of the organism [12-17], 
and comparative genomic analysis may contribute to systems synthetic biology or systems metabolic 
biology for targeting and engineering genetic circuits to create desirable cellular phenotypes. 
Transcriptome profiling uses DNA microarrays to decipher the expression levels of thousands of genes 
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under various biological conditions [14,15]. The results can be used to select candidate genes for 
modification based on systematic analysis of regulatory genes in response to genetic variations and 
environmental changes, or to identify novel factors for the enhancement of heterologous product 
secretion in metabolic pathways [64-68]. Proteome profiling is also useful in obtaining 
transcriptome-profiling data at the protein level. The metabolome comprises the entirety of information 
on metabolites present within and/or outside the cell under specified conditions [61-63]. It is expected to 
contribute significantly to the understanding of the cell and of synthetic circuit engineering in its 
metabolic pathways. In this paper, recent advances in the application of bioinformatics to systems 
synthetic biology and systems metabolic engineering through systems biology are reviewed using 
specific examples. 

Figure 1. The role of systems biology as an integrated platform in modern biological 
research Systems biology integrates information on genetics, proteins, DNA-protein 
binding, and metabolism with system dynamics modeling and system identification 
technology to develop models and mechanisms for the interpretation of phenotypes or 
behaviors of cellular physiology. Since large-scale data sets need to be mined, powerful 
computational tools are necessary. Based on system models and mechanisms in systems 
biology, synthetic genetic circuits are designed to investigate specific desired cellular 
behaviors of cellular physiology. Discrepancies between real and desired cellular behaviors 
are used as feedback to adjust system models and mechanisms. Systems biology is thus 
positioned to play the role of integrated platform for bioinformatics, systems synthetic 
biology, and systems metabolic engineering. 
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Although bioinformatic information (e.g., data on genetics, protein binding, and metabolism) is 
available, several stages of systems biology are required to help us understand via system dynamics 
modeling the underlying molecular mechanisms of genetic regulatory (GR) networks [18-24], 
protein-protein interaction (PPI) networks [25-29], and metabolic networks [61-63] under various 
biological conditions. At the first stage, a putative GR or PPI network is created by large-scale 
integration of knowledge such as information from publications and databases, and high-throughput data 
(from data mining or deep curation). Based on this network and dynamic modeling, the actual GR or PPI 
network of cellular physiology can be identified with system identification methods 
(reverse-engineering scheme [34]) by using specific microarray gene or protein expression data [ 1 8-28] . 
For example, GR networks (GRNs) have been constructed by dynamic modeling via microarray data for 
cell cycles [18,23,24], environmental stressors [28,44], photosynthesis [69], aging [34], and cancer [39]. 
PPI networks have been constructed for cancer [3,9,33,39,70], inflammation [41], biofilm formation [43], 
and infection by Candida albicans. Comparison of PPI networks between healthy and cancer cells can 
provide network markers for the investigation of the systematic mechanism of cancer [70]. The 
integration of cellular networks of GRs and PPIs provides deeper insight into actual biological networks 
and is more predictive than an approach without integration [71]. A systematic and efficient method to 
integrate different kinds of omics data for the construction of integrated cellular networks via microarray 
data have been provided based on coupling dynamic models and statistical assessments [44]. This 
method has been shown to be powerful and flexible, and can be used to construct integrated networks at 
different cellular levels to investigate cellular machinery under different biological conditions and for 
different species. Coupling dynamic models of the whole integrated cellular network is very useful for 
theoretical analyses and for further experiments in the field of network biology and synthetic biology. 

In short, synthetic biology is the engineering of biological systems to fulfill a particular purpose. It is 
possible to build living machines from off-the-shelf genetic devices by employing many of the same 
strategies that electrical engineers use to manufacture computer chips [47-49]. The main goal of this 
nascent field is the design and construction of biological systems with desired behaviors [51]. Synthetic 
biology envisions the redesign of natural biological systems as well as the construction of functional 
"genetic circuits" by using a set of powerful biotechniques for the automated synthesis of DNA 
molecules and their assembly into genes and microbial genomes [47]. Synthetic biology is predicted to 
have important applications in biotechnology, metabolic engineering, and medicine. It may 
revolutionize how we conceptualize and approach the engineering of biological systems [49]. As 
illustrated in Figure 1, synthetic genetic circuits can furthermore be employed to confirm network 
mechanisms derived using systems biology methods, and can be used as feedback for their improvement 
or revision. Synthetic biology is therefore expected to contribute significantly to a better understanding 
of the functioning of complex biological systems such as metabolic pathways. 

However, the development of synthetic gene networks is still difficult. Most newly created gene 
networks are nonfunctional because of intrinsic parameter fluctuations, environmental disturbances, and 
functional variations in the intra- and extracellular context. For this reason, the design method based on 
dynamic models for robust synthetic gene networks has become an important topic in synthetic 
biology [52-60,68]. These system-dynamics-based design methods for synthetic biology lead to systems 
synthetic biology. 
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Heterologous genes have previously been combined into pathways in metabolic engineering, 
generating a myriad of non-native biochemical products, including isoprenoids, hydroxyacids, biofuels, 
polypeptides, and biopolymers [64-67]. Synthetic biologists developed synthetic tools to engineer 
genetic devices capable of performing complex biological functions such as sensing cell states, counting 
cellular events, and implementing computational logic. These tools have been applied to the 
modification and control of metabolic pathways in several organisms. They consist of one or more parts 
that have been combined to perform a complex function, and provide metabolic engineers with novel 
ways of exerting cellular control over heterologous production pathways. Some synthetic biological 
devices with potential relevance in metabolic engineering include orthogonal inducible promoters, 
light-sensitive promoters, state sensors, spatiotemporal controllers, and logic gates, as well as promoter 
and ribosome binding site (RBS) libraries [36,37]. Since metabolic engineering seeks to control cellular 
metabolism and manipulate through heterologous pathways to maximize production of a desired 
molecule, metabolic engineers are need elegant methods for gathering bio-information about cells, their 
environment, and modulating gene expression in responses [37,61-63]. Hence, devices of synthetic 
biology promise to be a useful addition to the metabolic engineering toolbox. Some synthetic devices 
have already been used to increase product titers. However, many remain largely untested in an industrial 
setting, and the complexity of biology makes their application a feat of engineering [36,37,72]. From the 
systems biology perspective, continuous work with these devices can help elucidate design rules or aid 
the development of system dynamics models that facilitate their integration into metabolic industrial 
processes [36,37,72], and thereby lead to the development of systems metabolic engineering. 

Several mathematical techniques based on systems biology have been developed to analyze the 
systematic properties of complex biological networks. For example, system sensitivity of a biological 
network in response to various parameter variations has been analyzed to determine the systematic 
properties that affect the robustness and fragility of a biological network. System sensitivity analysis not 
only can reveal the robust stability of a biological network against various perturbations, but may also 
provide information about the controllability of a biological network [7-9]. The system response ability 
of a biological network is a measure of response to environmental signals or disturbances [28,34]. From 
the system theory perspective, robustness to intrinsic system variation and the ability to respond to 
external stimuli are two important and complementary system characteristics in the evaluation of system 
performance [33,34]. A more biological system that is robust toward a large amount of intrinsic 
parameter fluctuations is less responsive to environmental disturbances, and vice versa. A systems 
biology investigation of the aging-related gene network via microarray data found that network 
robustness increases and network response ability decreases during the aging process [34]. The 
sensitivity of a biological genetic system to environmental molecular noise is considered as an indication 
of the noise-filtering ability of the gene network [42]. Systems biology allows the measurement of its system 
characteristics, as well as the capabilities of the signal transduction pathway [8] Similarly, flux 
amplification of the metabolic pathway can be estimated, by using both system dynamics models. 

As nonlinear biological networks operate under different conditions of cellular homeostasis and 
homeodynamics, systems biology studies on complex biological models in the landscape of phenotypes 
are highly informative. These studies help discover possible equilibrium points (phenotypes) and 
dynamic behaviors, such as bifurcation, oscillation, robust stability, and phase drift to other equilibrium 
points (phenotype transition). Bifurcation analysis and phase-plane analysis of nonlinear dynamic 
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networks can be useful in predicting system behavior of biological networks under intrinsic parameter 
changes. Through systems biology approach and dynamic modeling [38-42], network robustness and 
noise filtering ability can be improved via feedback, redundancy, and modular schemes. This is why 
there are so many feedback loops, redundant genes, and modular structures at different scales of 
biological networks. A unifying mathematical framework based on nonlinear stochastic dynamic 
models [73-75] was recently proposed to describe different levels of stochastic biological networks 
under different parameter fluctuations, genetic variations, and environmental disturbances [29,30,59]. 
The phenotype robustness criteria of biological networks in systems, evolutionary, ecological, and 
synthetic biology were investigated from a systematic perspective on the basis of robust stabilization and 
filtering ability. Network robustness of biological networks can confer intrinsic robustness toward 
intrinsic parameter fluctuations, genetic robustness for buffering genetic variations, and environmental 
robustness for resisting environmental disturbances. It was found that if the sum of intrinsic robustness, 
genetic robustness, and environmental robustness is less than or equal to the network robustness, then 
the phenotype is robust in different levels of biological networks in systems, evolutionary, ecological, 
synthetic, and metabolic biology. These phenotype criteria at different levels of biological networks are 
useful for the design of synthetic and metabolic systems. A systems biology approach based on dynamic 
models can clearly provide not only a systematic insight into behaviors at different levels of biological 
networks, but also a design platform to improve system robustness, filtering ability, and transduction 
ability of synthetic and metabolic system networks, which are discussed in detail in the following sections. 

2. Systems Biology Approach to GRNs and PPI Networks via Bioinformatics Methodology 

2.1. Construction of GRNs via Microarray Data 

The construction of a GRN using a systems biology approach can be divided into two steps. The first 
step is construction of a potential (or putative) GRN from two high-throughput data sets, namely, a Chip 
network and a mutant network. The Chip data are based on experimental data published by Harbison 
et al. [14], who also used a genomic tiling array to identify the genomic region bound by transcription 
factors (TFs). The mutant data are the gene expression data matrix published by Hughes et al. [76] with 
different gene deletion mutants. In general, the potential GRN is suitable for all possible biological 
conditions. Therefore, the GRN for a specific biological condition needs to be confirmed using 
microarray gene expression data of the specific biological condition; that is, the real GRN is derived by 
pruning the potential GRN with specific microarray data. Let Xi(t) denote the gene expression of the i-th 
target gene in the potential GRN. The dynamic equation for the regulation of the i-th gene is then 
modeled as [21,22] 



where Xi(i) represents the mRNA expression level of target gene / at time t, Xj(i) represents the regulation 
function of the y'-th TF binding to the target gene /, ay denotes the regulatory ability from the y'-th 
regulatory gene to the z-th target gene (the positive sign indicates activation and the negative sign 
indicates repression), X t indicates the degradation effect of the present time point on the next time point 
t+l 9 ki represents the basal level, and «,•(*) is the stochastic noise due to model uncertainty and fluctuation 
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of microarray data of the target gene. Expression of the z'-th gene in (2.1) can be represented by the 
following regression equation [77] 



x,(r + l) = fo(0 ••• x N (t) !]• 



= #(03 +",(0 

where ^- denotes the regression vector, which can be obtained from microarray data. # z is the 
regulatory parameter vector of target gene /, which is to be estimated. 

The regulatory parameter can be estimated from microarray data of the corresponding gene through 
the least-squares or maximum-likelihood parameter estimation [77] (known as the reverse-engineering 
method). After regulator parameters in 6 t are estimated, the system order (i.e., the number of regulatory 
genes) is determined by model comparison (e.g., Akaike's information criterion (AIC)), and the P- value 
statistical method is employed to determine the significant regulatory genes for target gene /. This is 
done by pruning false-positive regulations in the potential GRN. That is, some ay is pruned because of 
false positive deletion. Based on the dynamic model in (2.1), the true GRN can then be constructed one 
target gene at a time through microarray data. Using similar methods, GRNs for yeast cell cycles [18,23,24], 
cancer cell cycles [78], stress response [44], and inflammation [41] can be constructed. 

2.2. Construction of PPI Networks 

The construction of PPI network via a systems biology approach is also a two-step process. The first 
step is constructing a potential PPI network via data mining from literature and databases such as 
BioGRID, SGD, and GO [16,17]. As this is only a candidate network based on many biological 
conditions, the second step is pruning false positive PPIs by using a protein expression microarray of a 
specific biological condition. For a target protein / in the potential PPI network, the dynamic model of 
protein activity is [19,20] 

M 

y t (t + 1) = y. (0 + 2 tyi «)yj (0 " PiVi (0 + h t + v > (0 (2.3) 

7=1 

where y££) represents the protein activity level of the target protein / at time t, by denotes the interaction 
ability of the y'-th interactive protein with the i-th protein, y/t) represents the activity level of the y'-th 
protein interacting with the target protein /, /?/ denotes the degradation effect of the protein, hi represents 
the basal activity level, and vfa) is the stochastic noise. The rate of PPI is proportional to the product of 
the concentrations of both proteins [7], that is, it is proportional to the probability of molecular collisions 
between them. The physical interpretation of Equation (2.3) is that the activity of target protein / at time 
t+\ is the combination of present protein activity, regulatory interactions with M interactive proteins, 
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levels of basal protein from other sources and M interactive proteins in the cell, and stochastic noise, less 
the protein degradation of the present state. The PPI dynamic equation of target protein i in the potential 
PPI network can be represented by the following regression equation [77] : 



y,y. 



M 



1] 



iM 

h 



+ v. 



(2.4) 



The interaction parameter 6 t can be estimated from protein profile microarray data by least-squares or 
maximum-likelihood parameter estimation [77] (if protein profile microarray data are unavailable, ten 
mRNA microarray data could be used to replace them, with some modification [19,20]). By using AIC 
to prune false positive interactions, the real PPI network can then be constructed one target protein at a 
time by following the above two-step procedure. Some dynamic metabolic pathways [20] and PPI 
networks of cancer [39] and inflammation [41] have recently been constructed by using the microarray 
data and AIC method. Comparison of PPI networks between healthy and cancer cells can provide 
network-based biomarkers for molecular investigation and diagnosis of cancer [70]. 

2.3. Construction of Integrated GRN and PPI Cellular Networks 

Living organisms have evolved complex mechanisms to respond to changes in environmental 
conditions. This is the case even in unicellular organisms like the yeast Saccharomyces cerevisiae. Such 
environmental changes, commonly termed as "stress", are harmful or even lethal to the survival of the 
cells, especially of microorganisms whose environment is very changeable. Cellular responses to 
stresses require complex cellular networks of sensing and signal transduction pathways, as well as 
metabolic pathways to adapt cell growth and proliferation for adjustments of gene expression programs, 
metabolic activities, and other features of the cell (see Figure 2). 

These regulatory systems and signal and metabolic pathways can therefore be suitably described by 
an integrated cellular network of transcription regulation and PPIs. The flowchart for constructing the 
integrated cellular network is shown in Figure SI. The same two steps already described apply 
(construction of a potential cellular network and pruning of false positives via specific microarrays of 
gene expression and protein expression data). Several kinds of omics data and database information 
were integrated, including data on microarray gene expression and protein expression, regulatory 
associations between TFs and genes, and PPI. From these data, the potential integrated GRN and PPI 
networks are retrieved. In the second step, based on the dynamic model of integrated gene/protein 
interactions in the potential cellular network, we find that [71] 
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y„ (t+i) = y n (0 + E ^ (0^ (0 + a n x„ (0 - (0 + a„ + v„ (0 



(2.5) 



where the regulation function z/t) is modeled as the sigmoid regulation function of yj(t), that is, for the 
protein activity profiles of transcription factor j [71], 

1 



z J (t) = f j (y j (t)) = 



l + exp{-(^(r)-u y )/t7 y } 



(2.6) 



where fjiyjif)) denotes the sigmoid function, uj and oj represent the mean and deviation of protein 
activity level of TF j 9 and a n denotes the translation effect from mRNAx^(0 to protein y n {f). 

Figure 2. Schematic diagram of the integrated cellular network. The integrated cellular 
network consists of two subnetworks. The signaling regulatory pathway (green) contains 
protein-protein interaction (PPIs). The gene regulatory network (yellow) contains transcription 
regulations. The transcription factors serve as the interface between the two subnetworks. 
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Interaction and coupling among genes and proteins in the integrated cellular network in (2.5) are 
described as follows. If some TFs y/t) at the end of the signal regulatory pathway regulate their target 
genes through the regulation function zj(t) = fjiyjif)) in (2.6), then the regulatory genes influence their 
corresponding proteins in signal and metabolic pathways through translation effect a n x n (t). The interplay 
between genes and proteins can be seen from their coupling dynamic equations in (2.5) and Figure 2. 
Here, the TFs serve as the interface between the signaling regulatory pathway and gene regulatory 
network. In other words, the interplay between transcription regulation and PPIs constitutes the 
framework of the integrated cellular network. 
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Based on the dynamic coupling equation in (2.5), the potential integrated cellular network can be 
linked through the regulatory parameter a$ between genes and their possible regulatory TFs and through 
the translation parameter a n for gene expression to protein expression. The potential signaling or 
metabolic pathways can be linked through the interaction parameter b nm between possible interaction 
proteins. Since omics data on the potential gene regulatory network and potential signaling or metabolic 
pathway only indicate possible TF-gene regulation and protein interactions, they should be confirmed 
using microarray data of gene and protein expressions. In particular, values of ay and b nm in (2.5) should 
be identified and validated by least-squares estimation via microarray data in a specific biological 
condition or phenotype [71]. Significant regulations and interactions between genes and proteins were 
detected using model selection methods such as AIC and statistical assessments like such as Student's 
t-test [77]. Based on the interface between gene regulatory and signal/metabolic networks (i.e., 
transcription factors) in a specific biological condition, the two networks are coupled together to form 
the integrated cellular network. The integrated cellular network for S. cerevisiae under hyperosmotic 
stress is shown in Figure 3. After the construction of GRN and PPI network from microarray data and 
bioinformatic method, various system characteristics of the biological network are estimated or 
measured using systems biology methods in the following sections so that these systematic 
methodologies can be applied the systematic design of systems synthetic biology and systems 
metabolic engineering. 

Figure 3. The S. cerevisiae integrated cellular network under hyperosmotic stress. By 
following the schematic diagram of an integrated genetic regulatory network (GRN) and 
PPI cellular network in Figure SI, the integrated cellular network of signaling regulatory 
pathway and GRN for hyperosmotic stress in S. cerevisiae is identified by dynamic 
modeling and data mining. Receptor proteins in the plasma membrane, signal regulatory 
pathways in the cytoplasm, and transcription factors and GRNs in the nucleus are used to 
construct an integrated cellular network for S. cerevisiae under hyperosmotic stress. 
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2.4. Network Robustness and Sensitivity Estimation via Microarray Data Using a Systems 
Biology Approach 

After the GRNs or PPI networks or their integrated cellular networks are constructed by dynamic 
modeling using microarray data, some characteristics of these biological networks need to be estimated 
to gain insight into their systematic mechanisms. "Robustness", the network ability to maintain systematic 
performance under intrinsic perturbations, and "response ability", the network sensitivity to respond to 
external stimuli or to transduce them to downstream regulators, are two important complementary or 
antagonistic system characteristics that must be considered when discussing biological network 
performance. However, these systematic features cannot be measured directly for all network components in 
the experimental procedure. Even biological processes such as development, differentiation, tumorigenesis, 
and aging are increasingly being described in terms of temporal ordering of highly orchestrated 
transcriptional programs [33,34]. The term robustness is encountered widely in diverse scientific fields, 
from engineering and control theory [73-75] to dynamic systems [79] and biology [80]. It is important 
to note that robustness describes a relative property, not an absolute one, as no system can maintain 
stability in all functions when it is perturbed. In other words, robustness is not immutable. 

Let the GRN of interest consist of N genes. After parameter identification via microarray data, the 
GRN in (2.1) can be represented by the following linear discrete-time dynamic system [33,34] 

x(t + 1) = Ax(t) + K + n(t) (2.7) 



A = 



a, 
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where the state vector x(i) = [x\(t) ••• xx(t)] T stands for the discrete-time mRNA expression levels of total 
N genes at times t = 1, 2, K. The system matrix A denotes gene interactions in the gene network 
estimated by microarray data, that is, dy denotes the estimated interaction of gene j with gene /. If / ^ y, 
then ki denotes the estimated basal level of the z-th gene. n(t) denotes the model residual and 
measurement noise. The steady state x s of x(t) is obtained as f— >oo: 



x -Ax 0 +K or x, 



K 



(2.8) 



To simplify analysis of "robustness" of the steady state (phenotype), the origin of the dynamic system 
is shifted to the steady state x s , that is, x{i) = x(f) + x s . This shift allows the following shifted dynamic 
system to be achieved by subtracting Equation (2.7) from Equation (2.8) [34]: 

x(t + l) = Ax(t) (2.9) 

Therefore, the robustness of the phenotype (i.e., the steady state) of the gene network becomes the 
robustness of the shifted gene network in Equation (2.9) at the origin x(t) = 0. In the following, its 
network robustness (tolerance to intrinsic network perturbation) is tested. Suppose the linear GRN 
suffers from intrinsic molecular perturbations mainly due to process noise, thermal fluctuation, or 
genetic mutations. The interactive matrix A is then perturbed as A(l + //), where r\ denotes the ratio of 
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intrinsic perturbation. That is, the corresponding additional system perturbation is AA = rjA [80]. A GRN 
with intrinsic perturbation can then be represented by 



Because quadratic stability with the Lyapunov energy-like function V(x) = x T {i)Px{i) > 0, with V(0) = 
0 for a positive symmetric matrix, P = P T > 0. The perturbative GRN in Equation (2.10) is robustly 
stable if AV(x) = V(x(t+X)) - V(x(t)) < 0, i.e., the energy of the gene network is not increased by intrinsic 
perturbations [79]. Based on this idea of robust stability, if the following inequality has a positive 
definite solution P = P T > 0 [34], 



then the perturbative gene network is robustly stable or the phenotype of the gene network is maintained 
under parametric perturbation ratio rj. The network robustness rf of the perturbative gene network is 
equal to the tolerance of the largest perturbation that does not violate network stability or still maintains 
the phenotype of the perturbative gene network: 

77° = max 77 

(2.12) 

subject to P>0, (2.11) 

That is, the network robustness is the maximum perturbation ratio rf tolerated by the GRN such that 
network stability (or phenotype) is still maintained. The constrained optimization problem in (2.12) can 
be solved by increasing rj until no positive solution P exists in Equation (2.1 1) up to the highest value 
possible without violating the robust stability in Equation (2.1 1). A positive definite solution P > 0 in 
Equation (2.1 1) can be easily obtained by using the linear matrix inequality (LMI) Toolbox of Matlab. 
This network robustness method has been used to measure the relative network robustness of multiple 
loops of a gene regulatory network associated with aging-related pathophysiological phenotypes by 
using previously reported microarray data. It profiles the effects of aging on gene expression in the 
thymus, spinal cord, and eye tissues in mice [34]. This aging-related GRN includes 16 genes (Figure 4). 
The relative robustness rf toward normalized perturbation for GRNs of young, aged, and 
calorie-restriction (CR) groups is shown in Table 1 . 

After the network robustness of aging-related GRNs is measured, the response abilities of genes to 
external stimuli are examined. Assuming that the GRN responds to external stimuli U(t), including 
upstream regulatory signals and external signals outside the network (e.g., hormones, carcinogens, 
oxidative stress, or ambient biomedical molecules) that could produce an effect on the gene network, the 
dynamic Equation (2.9) could be modified as follows: 



where U(t) = [u\(f) ••• uj^t)] T Represents external stimuli and Y(t) denotes the output signal response of 
specific genes of interest. For example, if the output signal response of the gene / to external stimuli U{t) 
is to be analyzed, then C= diag(0, 0, 0, 1, 0, 0, 0), i.e., all elements of C are zero except for a single 
element at the z-th diagonal component. For example, C = diag(0, 0, 1, 0, 0) describes the gene 




(2.10) 




(2.11) 



x(t + i) = Ax(t) + U(t) 
7(0 = Cx(t) 



(2.13) 
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response of the P53 gene. If the network response ability of the entire GRN to external stimuli is 
analyzed, then C = I (identity matrix). 

Figure 4. Multiple loops of a gene regulatory network associated with aging-related 
pathophyiological phenotypes. This aging-related gene regulatory network includes 16 
genes: FOXO s , NF-kB, P53, SIRT1, HIC1, Mdm2, Arfl, PTEN, P13K, Akt, JNK, IKK, IkB, 
BTG3, E3F1, and ATM. Dashed red lines and black arrows indicate negative and positive 
parameters of regulated interaction, respectively. 




The effect of input signals U(t) on output signal Y{t) is less than or equal to a positive value S, if the 
following inequality holds [79] 



jY T (t)Y{t) 



Yu T {t)U{t) 



<s 2 



(2.14) 



d denotes the upper bound of the effect of U{t) on Y{t) for all bounded input signals U{t). 8°, the smallest 
upper bound of S in (2.14), is called the "gene response ability (or sensitivity)" of the GRN. S° permits us 
to obtain a more systematic insight into the ability of gene response to external stimuli for individual 
genes or the entire GRN. From the system gain result in [79], the network response of the dynamic GRN 



has an upper bound 8 in (2.14), if there exists a positive definite P 
LMI [34]: 

" A T PA + C T C-P A T P 
PA T P-S 2 I 



<0 



P T > 0 solution to the following 



(2.15) 



That is, if the above LMI holds for some P > 0, then the effect of U(t) on Y(t) must be less than or 
equal to d, and Equation (2.14) holds. The response ability (minimal 5) of the GRN to external stimuli 
can be obtained by solving the following constrained optimization problem: 
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S° = min 5 

p>0 (2.16) 
subject to (2.15) 

This can be solved by decreasing the upper bound 8 in Equation (2.14) until no P > 0 exists in 
Equation (2.15) by using Matlab LMI Toolbox. By following Equation (2.16), the network response 
abilities of aging-related GRNs in Figure 4 can be measured and compared at different life stages. 
Results are shown in Table 1 . 

Table 1 shows that the "aged" group has higher network robustness but lower network response 
ability. In order to tolerate the intrinsic parameter fluctuations accumulated due to genetic mutation, the 
elderly GRN becomes more robust than the young GRN. Robust GRNs are less responsive to external 
stimuli; consequently, the network protective ability against external stimuli decreases in the older gene 
networks. However, a more robust GRN may harbor more accumulated genetic mutations, which 
through random drift might provide more evolutionary paths to other phenotypes of gene network and 
thus lead to some aging-related chronic diseases like cancers, metabolic disorders, and arthritis. Table 1 
also shows that the "young" network groups are less robust, with greater response to external stimuli. 
These GRNs are therefore also less robust to intrinsic perturbations and elicit a strong response toward 
external stimuli. This might imply that some gene expressions such networks could be easily 
reprogrammed to mediate downstream genes or regulators for further reactions. It also allows for 
modulation of gene expression in response to external stimuli, such as exposure to oxidative stress, 
carcinogens, and pro-inflammation molecules. These observations suggest that the gene regulatory 
network at the earlier stage of life or under conditions of CR may opportunely have protective 
adaptations to maintain intact regulatory structures and homeostasis of cellular functions. The purpose 
of the above analysis is to gain a greater understanding of the systematic protective and/or defensive 
mechanisms inherent to aging-related GRNs. The proposed robustness measurement methods may be 
used for future studies of GRNs involved in various biological processes and may have potential 
applications in gene therapy and drug target selection. It was also found that the network robustness of 
cancer cells is higher than that of normal cells and that the reverse is true for network response 
ability [33,34]. 

Table 1. The network robustness (rj°) and network response ability (S°) of a gene regulatory 
network with 16 aging-related genes (Figure 4) across different tissues at young, aged, and 
calorie-restrictive (CR) stages. 



Tissue 


Young 


Aged 


CR 


Thymas 








tf 


0.2310 


0.3750 


0.2050 


s° 


1.1653 


0.9463 


1.2279 


Spinal Cord 








rf 


0.2410 


0.6270 


0.1400 


S° 


1.1630 


0.9367 


1.2645 


Eye 








rf 


0.1600 


0.2910 




6° 


1.2376 


0.8798 
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2.5. Network-Based Biomarker Determination via Sample Microarray Data Using a Systems 
Biology Approach 

Biomarkers are used either in diagnostic evaluation to determine the health of a patient with or 
without a disease, or as a prognostic indicator to determine a patient's prognosis. Present biomarker 
identification methods, which strictly use gene expression profiles, cannot show how the different genes 
within the biomarker gene set are related to each other; that is, biomarkers are not identified from a 
systems biology perspective. Furthermore, the gene lists obtained for similarly diagnosed patients by 
different research groups differ widely and share few common genes. Here, a systems biology approach 
is introduced for the integration of microarray data and PPI information to develop a network-based 
biomarker for systematic investigation into the network mechanism of lung carcinogenesis and the 
diagnosis of lung cancer [70]. The network-based biomarker consists of two protein association 
networks constructed from cancer and noncancer samples. The proposed method can be widely applied 
to determining network-based biomarkers for other diseases. The overall flowchart of the proposed 
network-based biomarker approach is shown in Figure S2 [71]. First, PPI information and microarray 
sample data for smokers with and without cancer are used to construct potential PPI networks for cancer 
and noncancer samples. Since data for both samples are limited, the number of proteins used in potential 
PPI network construction is also restricted: to avoid overfitting in network construction, the maximum 
degree of proteins in the potential network are less than the number of samples, imposing a limit on the 
size of the potential network. 

Using a simple regression model, the potential PPI network is then further validated by the sample 
microarray data to highlight independent protein associations of both samples relative to their respective 
data sets. For a target protein / in the potential PPI network, the protein is described using the following 
protein association model [71]: 



(2.17) 



where y t {n) represents the gene expression level of the target protein / for the sample n, and a& denotes 
the association ability between the target protein y t {n) and protein yk{n) for sample n. N t represents the 
number of proteins interacting with the target protein /; it can be obtained from the rough PPI network. 
Si(n) denotes stochastic noise associated with other factors or model uncertainty. Equation (2.17) states 
that, biologically, the expression level of the target protein / is associated with the expression levels of 
interacting proteins. 

The associated interaction parameter a& in Equation (2. 17) is identified through maximum likelihood 
estimation [77] on microarray data. AIC and Student's t-test were employed for model order selection 
and for tests on the statistical significance of protein associations. Based on a#, two matrices are 
established to represent the cancer protein association network (CP AN) and the non-cancer protein 
association network (NPAN) as follows [77] 



C = 



&21,C ^22,C 



2K,C 



K\,C 



K2,C 



KK,C 



, N = 



a, 



12,N 



22, N 



'2a:,n 



a, 



X2,N 



KK,N 



(2.18) 
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where %c and %at indicate the quantitative protein association ability between protein / and protein j for 
CP AN and NPAN, respectively, and K is the number of proteins in the protein association network. The 
resulting CPAN and NPAN constitute the network-based biomarker used for identifying significant 
proteins in lung carcinogenesis through the diagnostic evaluation 

Y C =CY C +E C , 7 N =7V.7 N +£ N (2.19) 

where Y c = \y\,c(n) — yK,c(n)] T , Y N = \y\, N (n) — yK,N(n)] T denotes the vectors of expression levels, and£ c 
and E N indicate the noise vectors in cancer and non-cancer cases, respectively. A matrix indicating the 
difference between two protein association networks is defined as C - TV, i.e., 



dy y d^ 

d ' 2 1 d - 



12 



d, 



22 



IK 



V 2K 



d 



KK . 



^11,C ^11,N 
^21,C _ ^21,N 



^12,C ^12,N 
^22,C _ ^22,N 



®2K,C ~ ^2^,N 



a, 



K\,C 



a. 



7n,N 



a, 



K2,C 



X2,N 



a 



KK,C 



-a 



KK,N 



(2.20) 



where dy denotes the difference in protein association ability between CPAN and NPAN among proteins 
/ and j. Using matrix D to represent the difference in network structure between CPAN and NPAN, a 
carcinogenesis relevance value (CRV) was derived to quantify the correlation of each protein significant 
to lung carcinogenesis. To identify significant proteins, two important issues are taken into 
consideration. First, the magnitude of the association ability ay denotes the significance of association of 
one protein to another. A higher absolute value of ay implies that the two proteins are more tightly 
associated. Second, if a protein plays a more crucial role in lung carcinogenesis, then the difference in 
association numbers linked to the protein for CPAN and NPAN would be larger. For example, if one 
protein shares a strong association with many proteins in CPAN, but has weaker associations (no 
protein) in NPAN, then the protein in question is more likely to be involved in lung carcinogenesis. As a 
result, CRV is determined based on the difference in protein association abilities using the following 
equation [70]: 

CR y=ZKI (2.21) 

7=1 

For the /-th protein in the network-based biomarker, the implication of Equation (2.21) is that the 
CRV quantifies the extent of protein associations that differentiate CPAN from NPAN. 

The above network-based biomarker approach is applied to the molecular investigation and diagnosis 
of lung cancer. The primary data set of GSE41 15 (79 smokers with lung cancer and 73 smokers without 
lung cancer; obtained from the GEO database, http://www.ncbi.nlm.nih.gov/geo/) was used for 
construction of the network-based biomarker. The CPAN and NPAN in Equation (2.18), which consist 
of 399 and 393 protein associations respectively, constitute the network-based biomarker of lung cancer 
(Figure 5). The difference between CPAN and NPAN is shown in Figure 6. The CPAN indicates that 40 
identified proteins play significant roles in lung carcinogenesis (Table SI). 
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Figure 5. The constructed network-based biomarker. (A) Cancer protein association 
network (CP AN) obtained from C in (2.18) by maximum likelihood estimation, Akaike's 
information criterion (AIC) selection, and Student's t-test. (B) Non-cancer protein 
association network (NPAN) obtained from N in Equation (2.18) using the same criteria. 




Figure 6. The difference between CPAN and NPAN obtained from Equation (2.20) for 
network-based biomarkers for lung cancer. The significance of proteins (indicated by circle 
size) to the network-based marker is dependent on their CRVs in Equation (2.21), which are 
listed in Table S 1 . 
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2.6. On the Network Robustness and Filtering Ability versus Molecular Noise in GRNs Using a 
Stochastic System Approach 

Molecular noise has been shown to play many roles in the cellular functions of GRNs, including 
noise-driven divergence of cell fates and population heterogeneity, noise-induced amplification of 
signals, generation of errors in DNA replication leading to mutation and evolution, and maintenance of 
the quantitative individuality of cells [50]. Other cellular processes influenced by noise include 
ion-channel gating, neural firing, developmental modules, cytoskeleton dynamics, and motors [50]. 
Phase variation in pathogenic bacteria, in which cells alternate randomly between expressing certain 
genes and silencing others, is thought to be a form of cultivated noise [25]. These molecular-level noisy 
phenomena are deeply rooted in the statistical mechanical behavior of so-called nanoscale chemical 
systems, where concentrations of reactive molecule species are extremely low and, consequently, 
fluctuations (noises) in the reaction rates are large [50]. Even though the molecular fluctuations leading 
to phase variation seem random in the individual, regulatory factors tune the variation to ensure mean 
levels of heterogeneity for the population, i.e., the random molecular noises can be shown to be filtered 
or attenuated by the GRNs [25]. 

Since cellular molecular events are subject to significant thermal fluctuations and noisy processes 
with transcriptional control, alternative splicing, translation, diffusion, and chemical modification 
reactions, gene expression is best viewed as a stochastic process. Many observations suggest that 
molecular events underlying cellular physiology are subject to random fluctuations; these observations 
have led to the proposal of a stochastic model for gene expression and cellular functions [25]. Noise 
filtering can be considered from a signal processing perspective [74]. From this perspective, a pathway 
is viewed as an analog filter in terms of its frequency response. In terms of signal processing, these 
cellular pathways function as biological filters that transduce molecular signals and filter molecular 
noise [25]. 

For the convenience of illustration, the following linear biochemical dynamics of n genes GRN in 
Figure 7 [25] is simply considered initially: 

dx(t) 



dt 



= Nx(t) 



(2.22) 



where the concentration vector x{t) and stoichiometric matrix TV are given by 



x(t) = 



x 2 (t) 



N = 



TV 



TV, 



TV 



nl 



TV, 



TV 



where jt,-(f) denotes the concentration of the i-th gene, and Ny denotes the interaction between genes j 
and /. 

Suppose the linear GRN suffers intrinsic molecular fluctuations so that the stoichiometric matrix TV is 
perturbed as TV + ATV, 
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AN = 



'AN u (t) ... AN ln (t) 
AN nl (t) AN m (t) 



n(t) = Mn(0 



where ATV^ denotes the random parametric fluctuation of JV#; M// denotes the deterministic part 
(amplitude) of fluctuation; and n(t) is white Gaussian noise with zero mean and unit variance, and 
denotes the stochastic part of fluctuation, i.e., the stochastic part of fluctuation is absorbed to n{i). 

Figure 7. The linear n genes GRN with interaction Ny, intrinsic fluctuation ANy, gene 
expression and extrinsic fluctuation v t {t). 




The GRN with random parametric fluctuation can then be represented by 

dx(t) 



dt 



■ = (N + AN)x(t) 
Nx(t) + Mx(t)n(t) 



(2.23) 



For the convenience of analysis by the stochastic process, the above stochastic GRN can be 
represented as follows [74,75]: 

dx(t) = Nx(t)dt + Mx(t)da>(t) (2.24) 

where dco(t) = n{f)dt. co(t) is called the stochastic Wiener process or Brownian motion [73]. 

The stochastic dynamic equations of GRNs are always nonlinear. In order to meet the nonlinear 
stochastic regulatory networks, Equation (2.23) should be generalized as the following Langevin 
equation [25] 

dx(t) = N(x)dt + M (x) dco(t) (2.25) 
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where N(x) denotes the nonlinear interaction equation of the nonlinear GRN, and M(x)dco(i) is due to 
nonlinear intrinsic random fluctuation. 

Based on stochastic Lyapunov theory, let V(x) > 0 with V(0) = 0 denote the Lyapunov power-like 
function, then the stochastic GRN in Equation (2.23) or (2.24) is stochastically stable ifE[dV(x(t))/dt] < 
0 [73]. With the choice of V(x) = x T (t)Px(t) for some positive definite matrix P, the following result is 
derived. 

Proposition 1 [25]: 

The linear GRN with stochastic perturbation in Equation (2.23) or (2.24) is stochastically stable if the 
following LMI 

PN + N T P + M T PM < 0 (2.26) 

has a symmetric positive definite solution P > 0, i.e., the phenotype of the linear stochastic GRN is 
maintained under intrinsic stochastic fluctuation. 

Remark 1: (i) In the intrinsic noise-free case, the stable condition in (2.26) is reduced to PN + N T P < 
0, i.e., the eigenvalues of system matrix TV should be on the left-hand side of the complex domain. 
Obviously, if the LMI in Equation (2.26) has a positive solution P > 0, then the eigenvalues of N should 
be located on the far left-hand side of the complex domain with more negative real values in order to 
overcome the additional term M T PM due to intrinsic random noise, (ii) If some eigenvalues of system 
interaction matrix N are near the jco axis, then intrinsic random molecular fluctuations across the jco axis 
perturb these modes more easily such that the linear GRN becomes unstable. The LMI in Equation 
(2.26) can be rearranged to 

M T PM <-(PN + N T P) 



intrinsic robustness network robustness 



(2.27) 



that is, -(PN + N T P) in Equation (2.27) can be taken as a measure of network robustness, and M T PM due 
to the random parametric fluctuation can be taken as a measure of intrinsic robustness. The physical 
interpretation of Equation (2.27) is that if the network robustness can confer enough intrinsic robustness 
for tolerating intrinsic random parameter fluctuation, then the phenotype of the GRN is maintained. 

For the nonlinear stochastic GRN in Equation (2.25), the following result is also derived. 

Proposition 2 [25]: 

The nonlinear stochastic GRN in Equation (2.25) is still stochastically stable if the following 
Hamilton-Jacobi inequality (HJI) has a positive Lyapunov solution V(x) > 0 with V(0) = 0 

^0)Y A 7"/ \ 1 „T, ^ 2 V(X) 



dx 



N(x) + -M 1 (jc)- — y^M(x) < 0 (2.28) 
2 dx 



That is, the phenotype of the nonlinear stochastic GRN in Equation (2.25) is maintained under intrinsic 
stochastic fluctuation. Similarly, the HJI in Equation (2.28) can be rearranged as the following 
phenotype robustness criterion: 



M (x) Y^M(x) < - — N(x) 



2 dx 



dx 



(2.29) 



intrinsic robustness 



network robustness 



That is, if the network robustness of the nonlinear stochastic GRN in Equation (2.25) can confer 
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enough intrinsic robustness to tolerate intrinsic stochastic fluctuation, then the phenotype of the GRN 
is maintained. 

After the robust stability of GRN is guaranteed under intrinsic biochemical stochastic fluctuation, the 
effect of environmental random molecular noises on the GRN may be discussed. If the linear GRN in 
Equation (2.24) also suffers from environmental molecular noises v(t) = [v\(t) ••• v n (t)] T outside the 
network (see Figure 7), then 

dx(f) = (Nx(t) + Hv(t)) dt + Mx(t)da*(t) 

Z(t) = Cx(t) (2 ' 30) 

where H is a coupling matrix denoting the influence of environmental molecular signals v(i) on the 
GRN. Z(i) denotes the concentration of specific genes of interest. For example, if we want to examine 
the effect of molecular noises of co(t) and v(t) on gene / (i.e., x t {t)) 9 then we let C = diag(0...010...0). 
That is, every element of C is zero except for the i-th element. To investigate the effect of molecular 
noises on the whole GRN, then C = I, the identity matrix. The positive value p in the following 
inequality is then called the effect of environmental noises (or signals) on Z(i) in the stochastic GRN in 
Equation (2.30) withx(O) = 0 

E\~ Z T (t)Z(t)dt 

— <p (2.31) 



Jo 



v (t)v(t)dt 



That is, p is the upper bound of the effect of all environmental molecular signals or noises v(i) on the 
GRN. It is called the response level of the GRN. 

From this is derived the following response level result of the GRN in Equation (2.30). 

Proposition 3 [25]: 

The response level p of the linear stochastic GRN in Equation (2.30) is guaranteed if the following 
matrix inequality has a positive solution P > 0 

PN+N t P+M t PM+C t C+\pHH t P<Q (2 32) 

P 



By Shur complement [79], the above inequality is equivalent to the following LMI: 

<0 (2.33) 



PN + N T P + C T C + M T PM PH 



(PHY -p 

That is, if the above LMI has a positive solution P > 0, then the linear stochastic GRN in Equation (2.30) 
is robustly stable under random intrinsic molecular fluctuation and has a response level p to 
environmental molecular noises. 

The optimal response level po of the linear stochastic GRN (Equation (2.30)) can be obtained by 
solving the following constrained optimization 

p 0 = min p 

(2 34) 

subject to p> 0,P > 0 and LMI in (2.33) 



Remark 2: (i) If po < 1, then environmental molecular noises v(t) are attenuated by the GRN and po is 
called the filtering ability of the GRN, i.e., the GRN is less sensitive to environmental noises. If po > 1, 
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then environmental molecular noises are amplified by the GRN, i.e., the GRN is more sensitive to 
environmental noises, (ii) Substituting p 0 for p in Equation (2.32) and rearranging, the following 
phenotype robustness criterion of the stochastic GRN is derived 



AfPM + C T C+^PHH T P<-[PN+N T P) 
Po 



intrinsic robustness 



(2.35) 



network robustness 



environmental robustness 



The phenotype robustness criterion in (2.35) for the stochastic GRN in (2.30) can be denoted as 
"intrinsic robustness + environmental robustness < network robustness". In other words, if the sum of 
intrinsic robustness and environmental robustness is less than network robustness, then the phenotype of 
the stochastic GRN remains robust under the influence of stochastic intrinsic fluctuation and 
environmental noises. In order to maintain the phenotype robustness criterion, GRNs need negative 
feedback loops to improve the network robustness on the right-hand side of Equation (2.35). Parallel 
loops and modular and redundant structures are also required to reduce the effect of intrinsic parameter 
variations on the GRN and to resist environmental disturbance, i.e., to provide more intrinsic 
robustness and environmental robustness on the left hand side of Equation (2.35). This is why 
feedback loops, parallel loops, and modular and redundant structures are abundant in GRNs as they 
contribute to phenotype robustness and favor natural selection, (iii) If the network robustness in 
Equation (2.35) is not large enough to confer intrinsic robustness and environmental robustness 
simultaneously to maintain the phenotype of the GRN, some negative feedback gene loops are 
implemented as follows: 

dx (t ) = (( N + K ) x(t) + Hv(t) ) dt + Mx{t)dco{t) 
Z(t) = Cx(t) 



(2.36) 



where K denotes specific negative feedback loops that improve the network robustness of the 
stochastic GRN. In this situation, the phenotype criterion in Equation (2.35) is modified as 



M T PM +C T C+^PHH T P<- 



intrinsic robustness 



P(N+K)+(N+K) T P 



(2.37) 



environmental robustness 



network robustness 



It can be seen from the phenotype robustness criterion in Equation (2.37) that an adequate K can 
improve the network robustness of the stochastic GRN such that it better tolerates intrinsic stochastic 
fluctuation and is less responsive to environmental noise (i.e. smaller p<j). 

Similarly, the nonlinear stochastic GRN under environmental molecular noises (Equation (2.25)) 
should be modified as follows: 

dx(t) = (N(x) + H(x)v(t)) dt + M{x)dco(t) 
Z(t) = Cx(t) 



(2.38) 



The phenotype robustness criterion in Equation (2.37) is then modified as follows [25]: 



X Cs CsX ~\~ 



i 



4/>, 



'o V 



dV(x) 
dx 



dx 2 dx 



M<- 



dV(x) 
dx 



N(x) 



(2.39) 



environmental robustness 



intrinsic robustness 



network robustness 
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In general, it is very difficult to solve the HJI in Equations (2.29) and (2.39). Global linearization 
techniques [79,81] or T-S fuzzy methods [82,83] can be employed to interpolate several local linearized 
systems to approximate the nonlinear stochastic system in Equation (2.38) as follows: 

L 

dx(t) = Yj{ h i ( z ) N i x + H.v(t))dt + Mpcdcoit) 

i=\ (2.40) 

Z(0 = Cx(t) 

where h{z), i = l,---, L denotes the interpolation functions of global linearization or fuzzy bases with 
h,{z) > 0, Y}=i hi{z) = 1, and dx = (N t x + HiV)dt + Mjxdco is the z'-th local linearized GRN. 

In this situation, the network robustness criterion in Equation (2.39) can be modified accordingly [25] 

^C+^PKH^P+M^PM. <-(PN.+N T P), i = \,2,--;L 

p 0 2 ' ' -^v— - J ' I (2.41) 

v y / local intrinsic local network robustness 

local environmental robustness robustness 

From the phenotype robustness criterion in Equation (2.41), it is seen that if local network robustness 
is large enough for every local linearized GRN to simultaneously provide enough local intrinsic 
robustness for tolerating local random parameter fluctuations and enough local environmental 
robustness for resisting local environmental disturbance, then the phenotype of nonlinear stochastic 
GRN is maintained in spite of intrinsic stochastic parameter fluctuations and environmental molecular 
noises. If the phenotype robustness criterion in Equation (2.41) cannot be maintained, then negative 
feedback loops can be engineered to improve network robustness. These feedback loops may have 
potential applications for some types of therapy and drug target selection. 

Systematic methodologies for the analysis of system characteristics of GRNs and PPINs detailed 
above are useful for the systematic designs of systems synthetic biology discussed in the following section. 

3. Systems Synthetic Biology 

Synthetic biology can be expected to have important applications in biotechnology and medicine, and 
to contribute significantly to a better understanding of the functioning of complex biological systems. 
Synthetic biology is concerned with the engineering of biological systems that fulfill a particular 
purpose. This is achieved by transformative innovation that makes it possible to build living machines 
from off-the-shelf chemical ingredients, employing many of the same strategies that electrical engineers 
use to make computer chips. The main goal of this nascent discipline is the design and construction of 
biological systems with desired behaviors [45-49]. Synthetic biology envisions the redesign of natural 
biological systems based on a set of powerful biomolecular techniques for the automated synthesis of 
DNA molecules and their assembly into genes and microbial genomes, for greater efficiency, as well as 
for the construction of functional "genetic circuits" and metabolic pathways for practical purposes [49]. 
The construction of GRNs has recently demonstrated the feasibility of synthetic biology [84-87]. The 
design of robust gene networks still presents a difficult challenge, and most newly designed gene 
networks cannot function properly. Such design failures are mainly due to both intrinsic perturbations 
such as gene expression noises, splicing, mutation, evolutionary changes, and environmental 
disturbances such as changing extracellular environments [52,53]. Designing robust synthetic gene 
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networks that can tolerate intrinsic parameter perturbations, attenuate extrinsic disturbances, and 
function properly in a host cell is therefore an important topic in synthetic biology [55]. 

3.1. Design of Specifications-Based Systems Synthetic Biology 

Analysis of the dynamic properties of gene networks has been previously implemented using 
sensitivity analysis, either by qualitative simulation of coarse-grained models or by extensive numerical 
simulations of nonlinear differential equation models or stochastic dynamic models [56]. For 
applications in systems synthetic biology, however, these approaches are not satisfactory since local 
sensitivity analysis can provide only a partial description of all possible behaviors of a nonlinear gene 
network. In particular, it cannot guarantee that a synthetic network behaves as expected for all uncertain 
initial conditions, external disturbances, and parameter variations in a given range. Recently, Kuepfer 
et al. developed an approach based on semidefmite programming for partitioning parameter spaces of 
polynomial differential equation models into "feasible" and "infeasible" regions [68]. In this approach, 
"feasible" simply refers to the existence of a desired steady state of the synthetic network. Another approach 
using robustness analysis and tuning of synthetic networks was proposed by Batt et al. to provide a 
means to assess the robustness of a synthetic gene network with respect to parameter variations [51]. 
This approach allows searching for parameter sets for which a given property is satisfied, using a 
publicly available tool called RoVerGeNe. Several gene circuit design methods have recently been 
introduced insert or delete specific circuits in an existing gene network to modify its structure toward 
improved robust stability and filtering ability [38]. Robust synthetic gene network design, however, is a 
topic in itself, necessitating the design of a completely synthetic network with enough robust stability 
toward parameter fluctuations and with enough noise filtering ability to resist external disturbances, 
which allows it to function properly in a host cell. 

A robust synthetic biological design incorporating molecular noises has been developed based on 
stochastic game theory [52]. However, the intrinsic parameter fluctuations of synthetic gene networks 
have not been considered in the design procedure. In a previous study [53] some system design 
specifications (engineering designs) had been provided by users so that the designer must engineer an 
artificial gene network to meet these design specifications. For the convenience of problem description, 
a simple design example is provided to give an overview of design problems for robust synthetic gene 
networks. A more general problem scenario is treated in preparation for this overview. First, consider the 
cross-inhibition network shown in Figure 8. This network is synthesized with two genes (a, b) that code 
for two repressor proteins (A, B). More specifically, protein A represses the expression of gene b, and, at 
higher concentration, the expression of its own gene. Protein degradation is not regulated. This synthetic 
system can be modeled by the following synthetic equations: 



The state variables x a and Xb denote the concentrations of proteins A and B, respectively, k and y are 
the kinetic parameter and decay rate, respectively, r is the regulation function capturing the regulator 
effect of a transcriptional protein on gene expression, and has a smooth sigmoidal form (e.g., Hill 
function) [88]. 



K=h r al( X b ) r a2( X a)-7a X a 



(3.1) 



x b= k b r b( x a)-m 



(3.2) 
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Figure 8. A simple two-gene cross-inhibition network. The network's regulation functions 
are given in Equations (3.1) and (3.2). 




The simple cross-inhibition network in (3.1) and (3.2) can be represented by the following 
stoichiometric matrix equation [7] 



K o - r „ o 
o k„ o - 7b 



r al( X b) r al( X a) 
h( X a) 
X a 



(3.3) 



However, the stoichiometric matrix in vivo suffers from intrinsic parameter perturbations because of 
gene expression noises, splicing, mutation, evolutionary change, etc., as in [52] 

K^k a +Ak a n a r a ^r a +&r a n a 

k b ^k b + Ak b n b y b ^y b + Ay b n b 



(3.4) 



where Ak, and Ay, denote the amplitudes of fluctuations of the stochastic parameters and decay rates, 
respectively; and is a random white noise with zero mean and unit variance. 

Suppose the synthetic gene network also suffers from environmental disturbances because of 
changing extracellular environments and interactions with the cellular context in its host cell. The 
stochastic gene network can be then represented as 



k a +Ak a n a 
0 



k b +Ak b n b 



0 



0 



K 

0 



+ 



-Ya 
0 



0 



r b( X a) 
X a 



Ak a 
0 



- A Y a 
0 



r al( X b>a2( X a) 
r b( X a) 
X a 
X b 

r a^ X bYal( X a) 

x„ 



(3.5) 



0 0 

Ak b -Ay b ^ 



r b( X a) 



n b + 
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where x = [x a Xbf and v = [v a Vbf denote the state vector and environmental disturbance of the synthetic 
gene network in the host cell, respectively. These intrinsic parameter fluctuations and environmental 
disturbances may cause the engineered synthetic gene network to be dysfunctional in the host cell. 

If a synthetic gene network consists of n genes, then the stochastic gene network of Equation (3.5) in 
the host cell can be extended according to the following n-gene network dynamics [52] 

m 

x = Nf(x) + Y J M i g i (x)n i +v (3.6) 

i=\ 

where the state vector x = [xyx n ] T denotes the concentrations of proteins in the synthetic gene network, 
N denotes the corresponding stoichiometric matrix of the n-gene network, M/; / = 1,..., m denotes 
fluctuation matrices associated with independent random noise sources / = 1,..., m; and the elements 
of Mi denote the standard deviations of the corresponding parameter fluctuations, v = [vr--v„] r denotes 
the vector of external disturbances. This stochastic system is used to mimic the realistic dynamic 
behavior of a synthetic gene network of n genes in the host cell. As the network is subject to intrinsic 
parameter fluctuations and environmental disturbances in the context of the host cell, a robust synthetic 
gene network should be designed with the ability to not only tolerate parameter fluctuations and 
attenuate external disturbances, but also to achieve desired steady-state behaviors. 

For convenience of analysis and design, the stochastic dynamic in Equation (3.6) of a synthetic gene 
network in vivo can be represented by the following Ito stochastic differential equation [73,74] 

m 

dx = ( Nf(x) + v) + (xyHY t (3.7) 

1=1 

where Wi{t) is a standard Wiener process with dW l {t) = n t dt. 

To ensure correct and efficient operation of the gene network, several systematic design 
specifications should be imposed on it from the systematic engineering point of view. 

(i) The kinetic parameters and the decay rates in stoichiometric matrix should be chosen from a 
biologically feasible range 

Ne[N t ,N 2 ] 

(ii) The network should tolerate the stochastic kinetic parameter and decay rate fluctuations with 
prescribed standard deviations in Mi in the state-dependent noise term 

m 

i=\ 

(iii) The following desired steady state should be achieved 

x—>x d as t—>°° 

(iv) The following prescribed attenuation level p of environmental disturbance should be achieved 
(i.e., the Hao filtering ability) 

E L ( x ~ x df Q( x ~ x d) dt 2 

— 7^Z *P 2 (3.8) 

E\ v T vdt 
Jo 
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for all bounded environmental disturbances v(t), where Q > 0 is a symmetric weighting matrix and p is a 
prescribed attenuation level less than 1. That is, the effect of environmental disturbance v on the 
regulation error x - xj should be less than the attenuation level p from the average energy perspective. 

In order to achieve the desired steady state Xd and for convenience of design, the origin of the 
nonlinear stochastic system in Equation (3.7) should be shifted to x<i. Stabilizing the shifted nonlinear 
stochastic system at the origin would then also achieve xj 9 simplifying the design procedure. Let x =x- 
Xd. The following shifted stochastic synthetic genetic system is then derived [79]: 



in 

dx = ( Nf(x +x d )+v)dt+ ^Mfa (x + x d )dW t 



(3.9) 



1=1 



That is, the origin x = 0 of the stochastic system in Equation (3.9) is at the desired steady state Xd of the 
original stochastic system in Equation (3.7). N G [N1JV2] is then specified to tolerate the stochastic 
parameter fluctuation £j=i Mjgi(x + Xd)dWi and efficiently attenuate the environmental disturbance v to 

the prescribed level 



x T Qxdt 

Jo 



v T vdt 



(3.10) 



From the stochastic synthetic gene network Equation (3.9) in vivo, we obtain the following result. 
Proposition 4 [53]: 

If design kinetic parameters and decay rates in N E [N1JV2] are chosen such that the following HJI has 
a positive solution V(x) > 0 



^dV(x) \ , ^„ 1 fdV(x)y fdV(x)^ 

Nf(x + x d ) + xQx + 



v 



dx 



J 



4p 2 



V 



dx 



J 



V 



dx 



J 



1 



d 2 V(x) 
dx 2 



(3.11) 



M.g.(x + x d )<0 



then (a) the stochastic gene network in Equation (3.9) can achieve both the robust stabilization 
necessary to tolerate the intrinsic stochastic parameter perturbations and the prescribed attenuation 
level p of environmental disturbances, i.e., design specifications (i), (ii), and (iv) are satisfied, (b) If the 
stochastic gene network is free of environmental disturbances (v(t) = 0), then x— >0 or x— in 
probability, i.e., design specification (iii) is achieved. 

It is generally very difficult to specify N G [N\,N2] to solve HJI in Equation (3. 1 1) for V(x) > 0 using a 
systematic method. If all global linearizations are bound by a polytope consisting of M vertices 



dx 

dg x (x + x d ) 
dx 

dx 



Co 



M 



'\M 



K GmM j 



Vx 



(3.12) 
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where Co denotes the convex hull of a polytope with M vertices defined in Equation (3.12), then through 
the global linearization method [79,8 1 ], the state trajectories x{t) of the shifted gene network in Equation 
(3.9) can be represented by the convex combination of M linearized gene networks as 

M f m \ 

dx = Y j a j (x) NFjM + Yj M i G ij^ dW i +vdt ( 3 - 13 ) 



7=1 



V 



where the interpolation function affi satisfies 0 < aj(x) < 1 and X i — dj(x) = 1 . The trajectory of the 
nonlinear stochastic synthetic gene network in Equation (3.9) could thus be represented by the 
interpolated synthetic gene network in Equation (3.13). This yields the following result. 
Proposition 5 [53]: 

Assuming that design kinetic parameters and decay rates in N E [N\, N2] are chosen such that the 
following MLMIs have a common symmetric positive definite solution P > 0 



PNF. + F T N T P + V G T M T PM G + Q 



i=l 



P ~p 2 I 



<0, j = l 9 ... 9 M (3.14) 



then there are two results: (a) The stochastic synthetic gene network in Equation (3.7) is robustly stable 
toward intrinsic stochastic parameter perturbations and achieves a prescribed attenuation level p of 
environmental disturbances. That is, design specifications (i), (ii) and (iv) are satisfied, (b) If the 
synthetic gene network is free of environmental disturbance (v(t) = 0), then the synthetic gene network 
in Equation (3.7) may asymptotically converge to the desired steady state. (An in silico design example 
is shown in supplementary example 1.) 

Remark 3: (i) Gene circuit design can now be implemented using recombination technology to insert 
or delete TF binding sites in the promoter region of a regulated gene with the aim of increasing or decreasing 
the value of the kinetic parameter K t {i.e., different levels of affinity) of the regulated gene [89]. By 
inserting strong or weak binding sites, large or small values can be obtained. For example, the binding 
site of Ki = 1 is 1 0 times larger than that of k\ = 0.1 at the promoter region of target gene /. Changes to the 
decay rate y z can be implemented by shortening the 3 f polyadenylate tail (referred to as deadenylation), 
which primarily triggers decapping, resulting in 5 f to 3 ? exonucleolysis. Alternatively, removal of the 3 ? 
polyadenylate tail can increase y z [53]. Therefore, by shortening or elongating the gene's 3 f 
polyadenylate tail, we can increase or decrease the decay rate y t of gene /. Directed evolution methods are 
also useful in changing the elasticity (kinetic property of /q) and in designing biochemical circuits [53]. From 
a systems biology perspective, these advances in implementation techniques of K t and yi enable 
engineering of synthetic gene networks in the near future. 

(ii) Because of the uncertain values of v and initial state x(0) of the synthetic gene network (3.9) in the 
host cell, the following minimax design has also been considered to reset the uncertainty of v(i) and x(0) 
in vivo [52]. v{t) and x(0) are considered players maximizing their effects on regulation error x{t) when 
the kinetic parameters are specified as another player minimizing the effects of v and x(0) on x{t). 

E[ f x Qxdx 

min max (3.15) 

ke[k u k 2 } x(0) p^Tr V 7 
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A robust synthetic gene network based on dynamic game theory and fuzzy interpolation of local 
linearized linear systems can thus be efficiently designed [52]. 

(iii) When the desired output of the synthetic gene network is not a steady state but a time-varying 
signal like a specific binary or periodic signal, and the network structure becomes complex, the above 
systematic analysis methods are difficult to apply. Because natural selection of traits suited for 
environmental change is an important evolutionary mechanism, kinetic parameters of synthetic gene 
networks can be tuned by genetic algorithms (GAs) or evolution algorithms (EAs) to optimally track 
desired biological functions [54,58]. Based on the evolutionary network algorithm, the kinetic 
parameters may be tuned to maximize the fitness to some desired phenotype selected by natural 
selection. Consider the optimal tracking design of a synthetic gene network by network evolution shown 
in Figure 9. Let the tracking error be defined as 

<t)=y d (t)-y{t) 

The design purpose is then to tune the design parameter k t using evolutionary (or genetic) algorithms 
such that the network can achieve the following optimal tracking 

^l n ui E C eT (Oe(t)dt (3.16) 
where T p denotes the present time. Let 

J(k) = E\ Tp e T (t)e(t)dt (3.17) 
Jo 

Let the fitness function F(k) be defined as 

m= W) <318) 

Adapting a parameter vector (chromosome) k t G [k l L ,k l u ] by EA or GA to minimize J(k) for the 
desired network behavior tracking of the synthetic gene network is equivalent to maximizing the fitness 
function in Equation (3.18) to meet the natural selection. A robust biological network design with a 
desired output behavior y d(f) is therefore equivalent to as solution to the following fitness maximization 
problem using an evolutionary network method [54,58] 



F(k*)= max F(k) nim 

V J ke[k L ,k u ] V J ( 3 ' 19 ) 



The evolutionary network algorithm (or evolutionary systems biology algorithm [58]) is employed to 
solve the above fitness maximization problem via genetic operators such as selection, crossover, and 
mutation. It mimics natural selection in an evolutionary process to tune the kinetic parameter network of 
the synthetic gene network to an environmental change. A GA-based algorithm with binary coding of 
the chromosome has been proposed for the design of robust synthetic genetic oscillators with prescribed 
amplitude, period, and phase [54]. The oscillator is intended to allow protein concentrations to track 
desired periodic reference signals under intrinsic and environmental noises. Based on an evolutionary 
systems biology algorithm that encodes each chromosome in a real-valued vector, the design 
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parameters of target gene circuits can evolve to specific values in order to robustly track a desired 
biological function in spite of such interferences [58]. 

Figure 9. Block diagram of the optimal tracking scheme for synthetic biological circuit 
design using an evolutionary systems biology algorithm. Based on a network algorithm 
mimicking natural selection in an evolutionary process, the design parameters k of a 
synthetic biological circuit are tuned to minimize the tracking error between the desired 
logic circuit and the stochastic synthetic biological circuit, and to achieve the desired 
behavior tracking. 



Desired logic circuit 



+ 



Tracking 
error 



o 



Stochastic synthetic biological 
circuit in equation (3.9) 



y 



Tune the kinetic parameters 



Evolutionary 
systems biology 
algorithm 



3.2. Robust Synthetic Gene Network Design via Library-Based Search Methods 

Over the past decade, synthetic biology has made significant progress in designing biological parts. 
Even synthetic gene networks may have been constructed using a variety of biological components to 
achieve desired behaviors. A limitation on the development of complex synthetic gene networks 
intended to track specific reference trajectories is the lack of an efficient method for selecting suitable 
biological parts from libraries; experimental data in promoter libraries cannot be directly used for 
selecting adequate promoter parts. Current promoter libraries therefore need to be redefined based on 
promoter activity. This would allow development of library-based search methods. 

A dynamic model can be used in the indirect evaluation of the activities of promoter parts to help 
redefine existing promoter libraries. In the TetR-regulated promoter library shown in Figure 10, x(c, i) 
andX(c, t) respectively denote the concentrations of the mRNA and protein of gene yegfp (which is used 
to measure protein expression). The dynamic model of the promoter-regulation gene part is constructed 
with [56,57] 
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x(c,t) = P TetR (c,r)-]3x(c 9 t) 

X(c,t) = ax(c,t) - y yEGFP X{c,t) ^ 3 ' 20 ^ 

where fi and y y EGFP denote the degradation rates of mRNA and protein yEGFP, respectively, and a 
denotes the translation rate. The promoter regulation function P Te tR{c,r), which is dependent on repressor 
activity y and the choice of promoter c, has the form 

c — c 

P TetR (c,r) = c r + s r =c r +(c s -c r ) H TetR (r),c = [c r ,cje Lib TetR (3 2 1 ) 

l + \ r l K TetR) 6 

where the promoter c has two promoter activities c r and c s (the minimum and maximum values of the 
promoter regulation function PjetAc, r), respectively) for the TetR-regulated promoter library Lib Te tR- 
That is, c = (c r , c s ) G Lib Te tR- K Te tR and n Te tR denote the TetR-DNA binding affinity and binding 
cooperativity of regulatory protein TetR and DNA, respectively. H Te tR(r) is a Hill function capturing the 
effect of a regulatory protein. 

Based on the estimated promoter activities c = (c r , c s ) via maximum and minimum values of the 
steady state of protein expression data, some promoter libraries can be redefined in such a way that they 
can be efficiently selected from the design of the synthetic gene network (Table 2) [56,57]. Since a 
synthetic gene network always consists of a set of promoter-regulation gene circuits (Figure 10), the 
design of complex synthetic gene network addresses how to select a set of promoters from the 
corresponding promoter libraries that have promoter activities adequate for achieving the design 
specifications. A well-known gene circuit topology, the simple toggle switch, is shown for illustration 
purposes in Figure 1 1 . The toggle switch has two distinct stable states and can be reversibly switched 
between them by changing the inducers ATc and IPTG. Let x\(c\ 9 t), X2(c2, t), and x^ic^ t) denote the 
concentrations of mRNAs tetR, lad, and j/eg/p, respectively; and letXi(ci, t), Xi(c2, t), and X^{c^ t) 
denote the concentrations of proteins TetR, LacI, and yEGFP, respectively. The dynamic model of the 
toggle switch gene network in Figure 1 1 is then modeled as 

*1 ( C l > 0 = PLad ( C P r i ) - A ( C l > 0 

X x (c l9 t) = ax x (c x ,t)-y x X x (c x ,t) 

*2 { C 2 5 0 = PTetR { C 2 ? T 2 ) _ P X 2 { C 2 ? 0 

X 2 (c 2 ,t) = ax 2 ic 2 ,t)-y 2 X 2 (c 2 ,t) 

(3.22) 

^3 ( C 3 •> 0 = PTetR ( C 3 ? ^3 ) _ fl X 3 { C 3 ? 0 

X 3 (c 3 ,t) = ax 3 (c 39 t) — y 3 X 3 (c 3 ,t) 
y(c,t) = X 3 (c 3 ,t) 

c = , c 2 , c 3 ] , c x g Lib LacI ,c 2 ,c 3 £ Lib TetR 

The promoter regulation functions piacAcu r i), PTetR{c2, r 2 ), and pretR^, r 3 ) are dependent on the 
selection of promoters a, C2, and C3 from the corresponding promoter libraries in Table 2. The output of 
interest y{c, t) is dependent of the selected promoter set c=[c\, c 2 , c 3 ] with adequate promoter activities 
from the corresponding promoter libraries in Table 2. This dynamic toggle switch gene network model 
consists of three interactive dynamic models of promoter-regulation gene parts, as shown in 
Equation (3.20). 
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Table 2. Redefined TetR- and Lacl-regulated promoter libraries. The redefined TetR- and 
Lacl-regulated promoter libraries (/. e. , Lib Te tR and Lib Lac i) comprise different promoters (i. e. , 
Tk and k = 0,. . .,20) with their corresponding activities of c s and c r obtained from previous 
libraries of experimental data. 



TetR-regulated promoter library (Lib Te tR) 
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Promoter activity 
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Figure 10. Single schematic diagram of the synthetic promoter-regulation gene circuit. The 
existing TetR-regulated promoter library contains the minimum and maximum values of 
fluorescence [y z - m / w? yimax] corresponding to with and without TetR (repressor) binding. Based 
on the promoter regulation function (3.21) and these values, the promoter library is redefined 
for the design of synthetic gene networks (Table 2). 
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The multiobjective design approach for the H 2 /H^ synthetic gene network based on promoter 
libraries selects an adequate promoter set c=[c\ 9 c 2 , C3] from corresponding promoter libraries such that 
the following two design objectives are achieved simultaneously [56]: 

(i) Hoo desired noise attenuation level p& 



E^(y{c,t)-y r {t)) T Q{y{c,t)-y r {t)Yt 



<p\ force Libj, 7=1,..., 



E\y (t)v(t)dt 

(ii) H 2 optimal reference tracking: 

mm (y{c,t)-y r {t)f Q(y{c,t)-y r {t))dt 



m 



(3.23) 



eel 
j=l,.. J .,m 



(3-24) 



By solving an optimization problem with two constraints [56], adequate promoters can be selected 
from the corresponding libraries in Table 2 to achieve the two design objectives in Equations (3.23) and 
(3.24). Based on the synthetic toggle switch (Figure 11) and dynamic model (equation (3.22)) with 
intrinsic parameter variation and external disturbance, the adequate promoter set c = [a, c 2 , C3] [L 9 ,T 2 , 
Ls] is selected from promoter libraries in Table 2 to achieve the multi-objective H 2 IH^ reference 
tracking specified in Equations (3.23) and (3.24). Simulation results with v(t) = 10 x [ny-n^ 1 are shown 
in Figure 12. 

Figure 11. Synthetic gene circuit topology: simple toggle switch. The regulatory protein 
TetR, which is induced by ATc, inhibits the transcription of lad by binding promoter c 2 . 
TetR also inhibits transcription of yegfp by binding promoter C3 to repress the expression of 
the fluorescent protein yEGFP. The protein LacI, which is induced by the inducer IPTG, 
inhibits the transcription of tetR by binding promoter c\. The gene circuit has two distinct 
stable states, and can reversibly switch between them by changing the inducers ATc and 
IPTG. If an adequate promoter set c = [ci, c 2 , C3] is selected from corresponding promoter 
libraries, then yEGFP can be used to track the desired behaviors generated by a reference 
model. In the reference model, c\ is selected from the Lacl-regulated promoter library, and c 2 
and C3 are selected from the TetR-regulated promoter library in Table 2 (i.e., c\ E Libiaci, 
c 2 ,c 3 G Lib Te tR). 

IPTG 



tetR 






C l G Lib Lad 


TetR 




1 

ATc 






Cells 2013, 2 



668 



Figure 12. Simulation of toggle switch. By solving the LMI-constrained optimization 
problem of the H 2 IH^ design objective Equations (3.23) and (3.24) for the synthetic gene 
network in Figure 1 1 through the library searching method, an adequate promoter set c = [c\, 
ci, C3] = [Z9, T2, Zg] is selected from the corresponding promoter libraries. The inducer ATc 
is added to the synthetic gene network at 80-160 hours to induce the gene network, and then 
the inducer IPTG is added at 160-240 hours. The output y(c, t) clearly produces a robust 
track with the desired reference output y r (f). 
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Remark 4: Collective rhythms of GRNs, especially the synchronization of dynamic cells mediated 
by intercellular communication, have become a subject of considerable interest to biologists and 
theoreticians [90]. Synchronization of a population of synthetic genetic oscillators is an important 
consideration in practical applications, because a population distributed over different host cells needs to 
exploit molecular phenomena in a simultaneous manner in order to function as a biological entity. 
However, this synchronization of synthetic gene networks in different host cells may be corrupted by 
intrinsic kinetic parameter fluctuations and extrinsic environmental molecular noise. Therefore, robust 
synchronization is an important design topic in nonlinear stochastic coupled synthetic genetic oscillators 
with intrinsic kinetic parametric fluctuations and extrinsic molecular noise. A systems biology approach 
indicates [59,60] that if the synchronization robustness criterion is greater than or equal to the sum of the 
intrinsic robustness and extrinsic robustness, then the stochastic coupled synthetic oscillators can be 
robustly synchronized in spite of intrinsic parameter fluctuation and extrinsic noise. If the criterion for 
synchronization robustness is violated, then an external control scheme can be designed to improve 
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robustness by adding inducers to the coupled synthetic genetic network. These robust synchronization 
criteria and control methods are useful for a population of coupled synthetic networks with emergent 
synchronization behavior, especially for multicellular engineered gene networks [60]. 

4. Systems Metabolic Engineering 

Systems metabolic engineering aims to amplify or delete specific genes in metabolic pathways to 
perform metabolic engineering within a systems biology framework. Regulatory gene networks, 
metabolic networks, and other cellular networks can thus be engineered in an integrated system manner. 
Systems biological analysis via large-scale genome-wide analyses and computational bioinformatic 
tools can allow the rapid evaluation of the global physiology of a cell with respect to various cellular 
regulations. These include transcriptional and translational regulation, as well as metabolic engineering 
distribution [37,61-63]. Results of this systems biology approach can be used to predict targets for a 
metabolic engineering approach within the host strain. Furthermore, the integration of high-throughput, 
large-scale, genome-wide analyses with in silico simulation results might provide additional information 
on cellular status at various hierarchical levels from genome to fluxome. Strain selection for strain 
improvement by systems metabolic engineering is divided into several phases [37]: (i) A base strain is 
allowed to develop, (ii) The base strain is further engineered via synthetic biological methods, based on 
results obtained from high-throughput genome-wide bioinformatic data and systems biology 
computational analyses, (iii) The performance of this preliminary production strain is evaluated in actual 
fermentation, (iv) The results are then fed back into further strain development until a superior strain 
showing the desired performance is derived. 

By combining the results of multiple genome-wide analyses and computational analyses, systems 
biology may allow us to reach an unprecedented level of understanding of cellular physiology and 
characteristics [72], which can subsequently be used in a systems synthetic biology framework to design 
systems metabolic engineering strategies, for example, for strain improvement. However, novel 
computational methods such as large-scale data mining, multidimensional data integration, and 
data-driven network inference and deep curation schemes still need to be developed and applied to the 
integration and interpretation of heterogeneous large-scale genome-wide bioinformatic data, which are 
closely interconnected by complex regulatory and metabolic pathways [37,91]. Several robust 
biochemical circuit designs have recently been proposed to improve the network robustness of metabolic 
pathways. The proposed design schemes provide a systems biology method with potential applications 
in systems synthetic circuit design in systems metabolic engineering and systems drug design. Broadly, 
a metabolic network is a collection of enzymatic reactions that process cellular and intercellular 
metabolites. In systems metabolic engineering, the rates of reactions or fluxes, which correspond 
directly to changes in concentrations of substrates, enzymes, factors, or products, are often measured. 
Such systematic changes in concentrations can be expressed in terms of dynamic differential equations. 
The following S-system representation has been used for the last three decades as an efficient model for 
describing a dynamic metabolic network [7,27] 
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n+m n+m 



7=1 7=1 



X, = J] Xf - p t n * * =V t - V_ t , i-l, ...,n 



(4.1) 



7=1 7=1 



X n -a n Y[x^-P n \[x h ; 



7=1 7=1 



where Xi,---^+ m are metabolites, such as substrates, enzymes, factors, or products of a biochemical 
network, in which X\,---JC n denote n dependent variables and X n +\s~JC n + m denote the independent 
variables. In a metabolic network, intermediate metabolites and products are dependent variables, 
whereas substrates and enzymes are independent variables. The rate of change in X u X t is equal to the 
difference between two terms, one for production or accumulation, and the other for degradation or 
clearance. Each term is the product of a positive rate constant a z or fi t and all dependent and independent 
variables that directly affect production or degradation, respectively. Each variable Xj is raised to the 
power of a kinetic parameter gy and hy, which represents an activating effect ofXj onXi when its value is 
positive, and an inhibitive effect when its value is negative. V t and V- t represent aggregate flux into and 

out of the Xi pool. 

Construction of the S-system representation of a metabolic network and estimation of its kinetic 
parameters from experimental data are described in [7] and the literature references therein. The 
nonlinear parameter estimation problem of S-systems has recently been solved efficiently by evolution 
optimization methods [92]. It is generally difficult to study the network robustness or sensitivity of a 
nonlinear system such as Equation (4.1). Fortunately, many important characteristics of an S-system at 
or close to the steady state can be analyzed by using simple algebraic methods. Since most metabolic 
networks in nature operate near the steady state, at which inputs and outputs are almost balanced, the 
following focuses on the network robustness of metabolic networks at the steady state [7,27]. 

4.1. Robust Biochemical Circuit Design in Metabolic Networks 

Consider the steady state of the metabolic network in Equation (4.1), 



n+m n+m 



(4.2) 



7=1 7=1 



Taking the logarithm on both sides of Equation (4.2) and making some rearrangements, 



n n+m 



2^-^)lnX 7 .=ln^.-ln^.-2^-^)lnX., i = l,2,...,n 



(4.3) 




Introducing new variables and coefficients, 



yj =\nX_ 



7' 




(4.4) 



the steady state of the metabolic network is obtained as 
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(4.5) 



where 









V 








~a n • 






Y D = 


y n _ 


, b = 












■• a nn_ 


•4 



' ' ' a \,n+m 



In Equation (4.5), A D denotes the system matrix of the interactions between dependent variables Y D , 
and Aj indicates the interactions between dependent variables Y D and independent variables Yj. In the 
nominal parameter case, it is assumed that the inverse ofA D exists so that Y D can be solved uniquely, i.e., 
the metabolic network results in only one steady state (phenotype). Therefore, the steady state of the 
biochemical system is given by 



Y^A-^b-A^) 



(4.6) 



Suppose that parameter variations due to mutation, thermal changes, or disease can alter the kinetic 
properties of the steady state of a metabolic network as follows [27] 

(A D +AA D )(Y D +AY D )=(b+Ab)-(A I +M I )(Y I +AF 7 ) 



(4.7) 



Parameter perturbations are defined as 



AA D = 



Ab = 



Aa u 



Act 



Aa ;i 



Aa 
Ah 



Aa 



Ab 



, AA 1 = 



Ag„-Mi ••• A gu~ A K 

; Ag v -My ; 

A gl,n +i ~ ^n+l ■ ■ • A gl,n +m ~ ^n+m 

'■ A g U n + j- A hn + j '■ 

A Sn,n+\ ~ A K,n+l ' " ' A Sn,n+m ~ A K,n+m 



where AA D denotes perturbations due to kinetic parameter variations, Ab denotes perturbations due to 
rate constant variations, and AAj denotes perturbations due to kinetic parameter variations between 
dependent and independent variables. AA D can influence the existence of the steady state of the 
metabolic network. 

From Equation (4.7), we derive 

A D (I + A~ l AA D ) (Y D + AY D ) = (b + Ab) - (A, + AA, )Y r (4.8) 



If the following robustness condition holds [27] 



Km>| 2 <i 



(4.9) 



then the singular values of / + A D ~ ] AA D are free of zero and the inverse (/ + A D ~ l AA D )~ l exists. 
Therefore, the steady state of the perturbative metabolic network in Equation (4.8) is uniquely solved: 
as 



Y D + AY D = (I + A~ l AA D )~ l A- 1 [(b + Ab) -(A, + AA^Y, ] 



(4.10) 
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The above analysis implies that if the robustness condition in Equation (4.9) holds, then the steady 
state of a metabolic system is preserved under parameter variations AA D , i.e., Y D + AY D in Equation 
(4.10) has a small difference AY D from the nominal in Equation (4.6) under small perturbation. 
However, if the condition in Equation (4.9) does not hold, then individual values of / + A D ~ l AA D may be 
zero, the inverse (/ + A D ~ l AA D )~ l may not exist, and the steady state Y D + AY D may cease to exist under 
the parameter perturbation AA D . As an example, consider the singular value decomposition: 



(4.11) 



i=\ 



where a t denotes the z'-th singular value and u t , VjER" denote the corresponding left and right singular 
vectors, respectively. Therefore, if a parameter variation is specified as follows: 



M D =-a i ufv i , i = l,2,... 



(4.12) 



then 



A D +AA D =U J 



cr, 



(4.13) 



Obviously, the inverse (A D ~ l AA D )~ l or (I+A D ~ l AA D )~ l does not exist under the parameter perturbation 
in Equation (4.12). 

Remark 5: The parameter perturbations in the direction of singular vectors like Equation (4.12) are 
the network's points of fragility. The robustness prevents this kind of parameter variation and guarantees 
the existence of the steady state of the metabolic networks. When unexpected perturbations like 
Equation (4.12) are encountered, a catastrophic failure of the network follows. Robust circuit design is a 
necessary fail-safe mechanism in such situation. For example, the trehalose pathway in yeast consists of 
only a few metabolites that form a substrate cycle. It is governed by a surprisingly complex control 
system that is composed of several inhibiting or activating signaling mechanisms [7]. 

Equivalently, the network robustness of the metabolic network in Equation (4.9) can be rewritten as a 
more intuitive criterion for network robustness: 



(4.14) 



That is, ArA D T is the upper bound of AA D AA D T without violation of robust stability at steady state. If 
the network robustness criterion in Equation (4.14) holds, then the steady state of the perturbative 
metabolic network still exists. 

If a metabolic network cannot tolerate perturbations, i.e., the network robustness criterion in Equation 
(4.14) is violated, then robust control via a biochemical circuit design is a suitable remedy. Based on 
robustness analysis, we develop a biochemical circuit design scheme for the robust control of metabolic 
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networks. Consider the robust control system design of the metabolic network (in Equation (4.1)) by 
specific biochemical feedback circuits in a more general metabolic network form [27] 

n+m n n+m n 

X =®Y[ X U X ~P]\ X U X t i = l >-> n (4-15) 

7=1 k=\ 7=1 k=\ 

where X(* denotes a new biochemical control circuit with Xk regulating the production of X t by the 

kinetic parameter^. X\ k denotes a new biochemical control circuit with Xk regulating the degradation 

ofXi by the kinetic parameter l ik . The choice of regulating objects, X k and X h and the specification of the 
kinetic parameters,^ and are designed according to the feasibility of biochemical circuit linkage to 
achieve a desired robustness to tolerate AA D within the prescribed range of kinetic parameter 
perturbations in a metabolic network. Since fa and represent the elasticities of the corresponding 
enzymes in the designed control circuits, the implementation of control circuits is heavily dependent on 
the elasticity specification of these enzymes. 

Consider the robust control system of the metabolic network in Equation (4.15). By using a similar 
procedure to Equations (4.2)-(4.7) at steady state, 

(A D +F+M D \Y D +£Y D )=(b+/Sb)-(Aj +M I \Y I + AF 7 ) (4. 16) 

where fy and In are the kinetic parameters in F of the biochemical control circuit to be specified in Equation 
(4.15). 

Suppose we can find some F such that the inverse of (A D + F + AA D ) exists. Equation (4.16) is 
equivalent then to 

(A D + F) (I + (A D + FT 1 AA D ) (Y D + A Y D ) = (b + Ab) - (A, + AA, )(7 7 + A7 7 ) (4.17) 

Similar to Equation (4.14), a robust design scheme for the controlled metabolic system in Equation 
(4.17) is given by 

||(^ +^)- 1 A^|| 2 <1 or AA D AA T D <{A D +F){A D +F) T (4.18) 

In this case of robust circuit design, the design purpose is to specify feedback circuits such that the 
structural stability robustness of the metabolic network is improved, thus enabling the network 
robustness criterion in Equation (4.18) to tolerate larger parameter perturbations AA D . The phenotype 
(i.e., the steady state of the controlled metabolic network in Equation (4.16)) is then given by 

Y D + AY D = (I + (A D + F)~ l AA D ) _1 (A D + F)~ l [(b + Ab) - {A, + A^)(7 7 + AY,)] (4.19) 

A simple example of network robustness analysis and circuit design is given in supplementary 
example 2. 

4.2. Multipurpose Circuit Control Design of Metabolic Networks 

The robustness design above focuses on the tolerance of kinetic parameter perturbations AA D . The 
effects of rate constant variations Ab and of environmental changes or upstream regulatory changes A 7/ 
on output variations AY D should also be considered in the circuit design of metabolic networks to 
guarantee robustness against both intrinsic parameter variations and extrinsic environmental 
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perturbations. The sensitivity of A Y D to Ab in the designed metabolic network of Equation (4.16) is 
given by [7] 

A7 n 



D _ 



Ab 



= {A D +F)- 1 



The sensitivity of AY D to AYj is given by 



AY, 



AY, 



D _ 



= -(A D +F)- l A I 



(4.20) 



(4.21) 



It is more appealing to design a robust metabolic network with desired sensitivities to variations in 
rate constants and environmental signals, i.e., 



AY n 



Ab 



< s, 



AY n 



AY, 



(4.22) 



where the upper bounds si and S2 are prescribed in advance by the biochemical circuit designer. From 
Equations (4.20)-(4.22), the equivalent sensitivity criteria for Equation (4.22) are obtained as [27] 

I-sf(A D +F)(A D +F) T <0 
A, At -s 2 2 (A D + F)(A D +F) T <0 



(4.23) 



A multipurpose circuit design approach is therefore aimed to construct a biochemical circuit F to 
simultaneously satisfy the design requirements of network robustness in Equation (4.18) and network 
sensitivity in Equation (4.23). For example, suppose the goal is to design a robust biochemical circuit fn 
in supplementary example 2 to tolerate kinetic parameter variation AA D and satisfy the network 
sensitivity in Equation (4.22) or (4.23) with prescribed sensitivities of II AYu/Abh < IIAro/A&l^nominai = si 
= UdX = 3.42 and \AY D IAYh < IA7 D /A7J 2 ,nominai = s 2 = Un'Ufo = 2.66. Then/22 should be specified 
to satisfy the robust circuit design and the following inequalities: 
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(4.24) 



simultaneously. With the help of Matlab, the range of fz2 necessary to tolerate AA D in and satisfy the 
desired network sensitivity in Equation (4.24) is found to be [-1, -0.081]. We choose^ = -0.407 as a 
design example, which is a negative self-regulation. It has been found to efficiently eliminate the effect 
of parameter variations by negative compensation. About 10% of yeast genes encoding regulators are 
negatively self-regulating; thus, this mechanism seems to be important for maintaining robustness in 
yeast [27]. The metabolic network and time responses are shown in Figure S5(d) and (e), respectively. 

Similarly, suppose the goal in the design case in equation is to specify fn and hi to satisfy equations 
and (4.24), i.e., 
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(4.25) 



where s\ and S2 are the same as above. Similarly, the ranges of fn and I22 necessary to tolerate AA D in 
equation () and satisfy the desired sensitivity criteria in Equation (4.25) are found to be [-1, 0] and 
[0,1], respectively. /12 = -0.08 and / 22 = 0.31 are chosen as a design example. The metabolic network 
and time responses are shown in Figure S5(f) and (g), respectively. Another example of a TCA cycle is 
given in supplementary example 3. 

Given that the goal is to design a robust metabolic circuit to tolerate the perturbation AA D and to 
achieve the desired sensitivities s\ and s 2 in Equation (4.23), the robust control circuit design problem 
can be reduced to specifying F to satisfy the following multipurpose control circuit design derived from 
Equations (4.18) and (4.23): 

AA D AA T D <(A D +F)(A D +F) T 

I- s f(A D +F)(A D +F) T <0 (4.26) 
A I Af-s 2 2 (A D +F)(A D +F) T <0 

with nominal sensitivities s\ = IA D ~ 1 \\ 2 = 8.3685 and S2 = lA D ~ l Ai\\ 2 = 7.5464. 

If one biochemical control pathway with kinetic parameter fn is designed to satisfy the 
multi-objective design criteria in Equation (4.26), then the range of fn is found to be within [-0.8, -0.1]. 
If fn is chosen as -0.2 (Figure S6(a), blue line), then the time responses of the designed TCA cycle 
network shown in Figure S6(d), which match the desired properties of the proposed design method, are 
obtained. That is, the robust controlled biochemical network not only can tolerate AA D (to preserve its 
phenotype under parameter perturbations) but also retains sensitivity to environmental molecules in the 
nominal case. If dynamic circuit design is employed to implement the biochemical circuit fn, then an 
enzyme capable of catalyzing the reaction X2~^X\ is required. Additionally, a TF (Z) has to be found 

such that oxaloacetate2 (X2) can bind to the promoter of the enzyme's inhibitor gene, as fn is negative. 
The concentration of X 2 could therefore regulate X\ through the kinetic parameter f 2 . The elasticity of 
the enzyme inhibitor's gene sequence has to be modulated then to the specified performance by rational 
design or directed evolution. This allows the construction of the biochemical control circuit fn- 

The previous section discusses robust circuit design of metabolic networks based on the steady state 
of an S-system model. In the following section, a stochastic dynamic model approach is taken based on 
the systems biology approach outlined in Section 3. 

4.3. Robust Control Circuit Design of Stochastic Dynamic Metabolic Networks 

Based on sampled data from biochemical experiments and the delayed effect of molecular diffusion 
and transport in cells, the linear metabolic regulatory network can be suitably modeled by the following 
discrete-time dynamic system: 
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X(t + l) = AX(t), X(0) = X 0 
y{t) = CX{t) 



(4.27) 



where the state vector X(t) denotes the discrete-time expression vector of the molecules (mRNAs, 
proteins, or other chemical complexes in the biochemical regulatory network) at time t, A denotes the 
stoichiometric (interactive) matrix among these molecules (see Figure 13), and y(t) denotes the 
metabolic output. Thus, 

x 2 (t) 



X(t) = 



x,(t) 



A = 



a, 



*nl 



where / = n are molecular concentrations of mRNAs, proteins, and other complexes in the 
biochemical network. 



Figure 13. Linear metabolic network of n molecules with intrinsic parameter fluctuation Aay 
and extrinsic noise Vify denotes the biochemical circuit design from Xj to x z to improve 
network robustness stability and noise-filtering ability. 




Cells 2013, 2 



677 



Suppose the kinetic parameters of a biochemical regulatory network of metabolic processes are 
affected by the following intrinsic perturbation and environmental disturbance: 



v i-=i y 

y(t) = CX(t) 



X{t) + B v v{t) 

(4.28) 



where Xf=i ^/«i(0 denotes the intrinsic parameter fluctuations due to an Z random fluctuation source 
(e.g., thermal fluctuation, alternative splicing, molecular diffusion, etc.), and rii(t) denotes the z-th 
random noise with the statistics E[w z (^)] = 0 and E[w z 2 (^)] = cr z 2 , / = l,---, L. v{i) denotes the 
environmental disturbance. 

Based on the robust stabilization and disturbance filtering described in Section 2, the environmental 
disturbance attenuation level p of the metabolic network is denoted as 

Ho| 2 



v 



(OIL 



< P (4.29) 



Following the systems biology approach in Section 2, the following robust result for a linear 
stochastic metabolic network is derived. 
Proposition 6 [28]: 

For a linear stochastic metabolic network with intrinsic parameter fluctuation and environmental 
disturbance, if there exists a positive definite matrix P = P T >0 such that the following matrix inequality 
holds for a desired disturbance attenuation level p 

L 

A T PA -P+Y, a ? A ? PA i + cTc + ( ATpB v f ( B v PB v - P 2l T l ATpB v < 0 (4.30) 

i=\ 

then the intrinsic parameter perturbation can be tolerated and environmental disturbance can be 
attenuated to a level p in the stochastic metabolic network Equation (4.28). 
The phenotype robustness criterion in Equation (4.30) could be rewritten as 

L 

Y j e- A -P A i + C T C + (A T PB v )\B T v PB v -p 2 iy x A T PB v < P-A T PA 



vill ^ , environmentdrobustness network robustness 

intrinsic robustness 



(4.31) 



Remark 6: (i) The physical interpretation of the phenotype robustness criterion of the metabolic 
network in Equation (4.31) is that if intrinsic robustness allowing tolerance of intrinsic parameter 
fluctuation and environmental robustness, as well as attenuation of environmental disturbance are 
simultaneously conferred by the network robustness, then the phenotype of the metabolic network is 
maintained. It can be shown that if the eigenvalues of A are closer to the origin, then the network 
robustness is commensurately larger. 

(ii) To solve the phenotype robustness criterion in Equation (4.30), it can be transformed into the 
following equivalent LMI: 



A T PA-P+C T C+Y a? A? PA i A T PB 

/ J III V 



B T V PA B'PB-p'I 



< 0 (4.32) 
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According to its definition, the disturbance attenuation (filtering) ability po (i.e., the minimum p) can 
be obtained by solving the following constrained optimization problem: 



Po =mm/7 



subject to P > 0, and LMI in (4.32) 



(4.33) 



If the metabolic network in Equation (4.28) cannot achieve a prescribed disturbance attenuation 
ability pp (i.e., pp < po) specified for therapeutic or biotechnological purposes, then a biochemical circuit 
design using state feedback is necessary: 



X(t + l) = 



X(t) + B v v(t) 



(4.34) 



A + F + ^A^t) 
y{t) = CX(t) 

The phenotype robustness criterion in Equation (4.31) is then modified as 

X afAfPA, + C T C + (A T PB v f (B T v PB v - p 2 p iy l A T PB v < P — (A + F) T P(A + F) (4 J5) 



environmental robustness 



network robustness 



intrinsic robustness 



The negative feedback circuit F could shift the eigenvalues of A to the origin, thereby improving the 
network robustness of the metabolic network on the right-hand side of Equation (4.35) and achieving the 
prescribed attenuation level pp < po. In general, the interactions of a biochemical regulatory network in 
metabolic processes are nonlinear in real biosystems. In this situation, a nonlinear biochemical 
regulatory network of metabolic pathways under intrinsic stochastic parameter perturbation and 
environmental disturbance can be represented based on the stochastic dynamic model of systems 
biology in Section 2: 

X(t + 1) = f(x) + f i f i (x)n t (t) + B v v(t) 



(4.36) 



y(t) = CX(t) 



In the following, the robust stability and filtering ability po on v{k) at y{k) in the stochastic nonlinear 
metabolic network in Equation (4.36) is discussed. 
Proposition? [28]: 

If the following matrix inequality holds for a Lyapunov function V{x{k)) > 0 and a noise attenuation 
level p 



dX(t) 



(f(X(t))-X(t)) + 



dX(t) 



dx\t) 



2P 1 



dX(t) 



B B' 

v v 



dX(t) 
dX(t) 



+\if i T (*(t)) 

^ 1=1 

+ -X T (t)C T CX(t)<0 



(4.37) 



then the effect of environmental disturbance v{k) on y{k) is less than p, i.e., the robust disturbance 
attenuation level p is achieved for the nonlinear metabolic network in Equation (4.36). 

The filtering ability po of the nonlinear stochastic metabolic network in Equation (4.36) that 
attenuates environmental disturbance can be obtained by solving the following constrained optimization 
problem: 
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p 0 = mm p 

V(X)>0 

subject to HJI in (4.37) 



(4.38) 



Following the robustness result discussed in Section 2, the phenotype robustness criterion for the 
nonlinear stochastic metabolic network in Equation (4.36) is given as 

12 T 



i^ r avwo) 



intrinsic robustness 



+- 



dX(t) 



dX(t) 



+- 



2P 2 



dX(t) 



B.B. 



3FXT(0) 
dX(t) 



+ -X T (t)C T CX(t) 



(4.39) 



environmental robustness 



<- 



'dV(X(t)) V 



dX(t) 



{f(X(t))-X(t)) 



network robustness 



Suppose the network robustness of the metabolic network cannot confer enough intrinsic robustness 
and environmental robustness to maintain the phenotype of the metabolic network. In this situation, 
some negative feedback loops Fg(X(i)) should be implemented to improve network robustness as 
follows: 



(4.40) 



X(t + l)= f(X(t)) + Fg(X(t)) + £ f t (X(t))n t (t) + B v v(t) 

i=\ 

y(t) = CX(t) 

In this case, the phenotype robustness criterion in Equation (4.39) should be modified as follows: 



i=\ 



dx\t) 



intrinsic robustness 
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+- 

2 



r dV(X(t))) T ( dV(X(t)) 



v 



dX(t) 



dX(t) 



+- 



2p 2 



dX(t) 



B.B T 



dX(t) 



environmental robustness 



< - 



r dV(X(t)) V 



dX(t) 



{f(X(t)) + Fg(X(t))-X(t)) 



network robustness 



That is, the negative feedback loops Fg{X) are employed to improve network robustness to tolerate 
more fluctuations of the intrinsic parameter and to filter more environmental disturbances. 

It is generally very difficult to solve the HJI in Equations (4.37), (4.39), and (4.41) for robust 
chemical circuit design of metabolic networks. Based on the global linearization method in Equation 
(3.12) [79,81] or fuzzy interpolation methods [83,93], the nonlinear stochastic metabolic network can be 
approximated by interpolating M local metabolic networks as follows: 
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(4.42) 




^(0 = CAT(0 



where f[X) 9 g(X) 9 and yp) are approximated by Y£\ a/X)AjX(t) 9 ^j£ x GjX(t), and ££i ^(0, 
respectively. In this situation, the following robust chemical circuit design for a metabolic network 
with a prescribed noise attenuation level pp is derived. 
Proposition 8: 

Suppose the negative feedback loop F is designed for the metabolic network in Equation (4.40) or 
(4.42), such that the following LMIs have a positive solution P > 0: 



Thus, the prescribed disturbance attenuation level is achieved by the biochemical circuit design. 

If the optimal robust filtering design is employed for the nonlinear stochastic metabolic network in 
Equation (4.40) or (4.42), then the control matrix F is specified to solve the following constrained 
optimization problem: 



where po in Equation (4.44) is the disturbance (noise)-filtering ability of the optimally controlled 
metabolic network in Equation (4.40) or (4.42) (see supplementary example 4). 

This robust biochemical circuit design example in a metabolic network illustrates that the proposed 
systematic method of biochemical circuit design produces metabolic networks that not only tolerate 
more fluctuations of the intrinsic random parameter but also filter more environmental disturbances. 
Systems metabolic engineering design is thus capable of improving the network robustness of metabolic 
networks and efficiently attenuating the effect of environmental noises. This approach may serve as a 
basis for drug designs against genetic perturbations, pathological environmental cues (such as infectious 
agents or chemical carcinogens), or both. 

Network robustness is a systematic property that allows a metabolic network to maintain its 
biochemical function or to generate biochemical products despite external disturbance and intrinsic 
parametric fluctuation. It is from a class of fundamental and ubiquitously observed systems-level 
phenomena that cannot be understood by observing individual components only. Thus, a study of network 
robustness at one level is simply a study of why, when, and how metabolic networks function properly or 
otherwise [47,80,94]. If network robustness is not large enough to obtain sufficient intrinsic robustness 
to tolerate random parameter fluctuations and environmental robustness to filter environmental 
molecular noises, then a robust biochemical circuit design is necessary to improve robustness such that 
the network can maintain its function or phenotype. The proposed robust circuit design principles could 
potentially be used for robust biosynthetic network design with applications in drug design, gene 
therapy, and metabolic engineering. Future focus may be on robust chemical circuit design for such 
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applications or construction of other pathways by nanotechnology and metabolic engineering [28]. If 
network robustness, intrinsic robustness, and environmental robustness are taken into consideration in 
metabolic network engineering, such networks would function more reliably and efficiently. 

With the exception of implementing feedback biochemical circuits in metabolic networks, systems 
synthetic biology techniques detailed in Section 3 could be combined with systems metabolic 
engineering approaches to engineer complete metabolic pathways or networks to produce biochemical 
products that could not be produced in the host cell. As an example, E. coli does not have 
non-fermentative pathways for isobutanol (see Figure 14(A)). Suppose a synthetic pathway mE. coli is 
required for the production of butanol as biofuel. In this application, it is necessary to engineer 
transcription and translation genes to produce the enzymes AlsS, IlvC, IlvD, and Kdc (KIVD), which are 
necessary to catalyze the isobutanol metabolic pathway (Figure 14(B)). The last enzyme of the 
metabolic pathway (Adh) is already present in E. coli. Promoters and RBSs with different regulatory 
abilities can be selected from the corresponding libraries. With adequate computation and simulation 
following the systems biology approach (see Section 2), a synthetic metabolic pathway could be 
designed and engineered (see Section 3) to achieve high-yield and high-specificity metabolic production 
of isobutanol. 

Figure 14. Engineered synthetic metabolic pathway for isobutanol production in E. coli. 
(A) Schematic representation of engineered isobutanol production pathway. (B) Engineered 
synthetic genetic circuit to generate the enzymes necessary for pathway in (A) for isobutanol 
production in E. coli. 
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5. Future Challenges in Systems Biology 

Systems biology approaches to synthetic biology and metabolic engineering are built on molecular 
and genetic findings, and results of studies in omics fields such as genomics, proteomics, and 
metabolomics. The main concerns faced by systems biology are the complexity and dynamic character 
of biological systems, the vast quantities of biological data, and the fragmented nature of biological 
knowledge. These fragments of information need to be integrated with nonlinear stochastic dynamic 
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models at different system levels by suitable computational tools before they can be applied to system 
synthetic biology and systems metabolic engineering. 

The first challenges of systems biology are how to enrich omics data (e.g., updating and improving 
protein microarray, ChlP-chip databases, and pathway database for systematic studies discussed in 
chapters 2-4) and develop more powerful computational tools for sophisticated data handling, advanced 
modeling, integrated analysis, and knowledge integration. A systems biology approach based on the 
integration of these enriched omics data and advanced computational models will enable us to predict 
the behaviors of biological systems more precisely. This will increase our understanding of the 
underlying molecular mechanisms and our ability to efficiently predict the effects of designed genetic 
circuits in metabolic engineering as well as the impact of perturbations on biological systems in 
drug treatment. 

The second challenge is how to increase the capability of researchers to navigate and relate various 
data and knowledge resources using the integrated platform of bioinformatics, systems synthetic 
biology, and system metabolic engineering to enable innovations. Connecting genomics, molecular 
networks, and physiology will provide us with deeper understanding of how individual differences in the 
genome affect physiological processes through alterations in molecular networks. Therefore, predictive 
and preventive medicine based on network-based biomarkers inevitably lead to personalized medicine 
that may revolutionize healthcare in the future [7]. The third challenge of systems biology is how to 
integrate bioinformatics, gene circuit design, and metabolic engineering technologies for systems drug 
design in future predictive and preventative medicine. 

Drastically increasing oil consumption, exhaustion of natural resources, and global warming are 
worldwide concerns at present. They have spurred research in the areas of bioenergy and biorefining, as 
well as microbial production of bulk chemicals and materials. Systems biology will play an important 
role in these undertakings, collaborating with systems synthetic biology and systems metabolic 
engineering to generate bio-based products to replace some, if not most, currently used chemicals and 
materials. The fourth challenge of systems biology is how to make these bio-based processes competitive. 
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